Main

We used imputation to assess association with T1D across 2.6 million polymorphic SNPs from the International HapMap Project in a total of 7,514 cases and 9,405 controls of European ancestry from three existing GWA studies: Wellcome Trust Case-Control Consortium (WTCCC; UK)2, Genetics of Kidneys in Diabetes–National Institute of Mental Health (GoKinD-NIMH; USA)3 and Type 1 Diabetes Genetic Consortium (T1DGC; UK)4 (Supplementary Table 1). We used the R package snpMatrix7 to conduct the imputation and calculate single SNP association score tests for each HapMap SNP. The score tests were based on the Cochran-Armitage test, with a Mantel extension to allow combination over different strata (UK region in the case of the WTCCC and T1DGC samples, and estimated ancestry score derived from principal components in the case of the GoKinD-NIMH samples3). For imputed SNPs, we calculated the score statistics using the expected value of the imputed SNP, given observed SNPs, with the expectation calculated under the null hypothesis.

Overall, there was some overdispersion of test statistics (λ = 1.14 and 1.09 for 1 and 2 degrees of freedom, respectively). This was consistent with the large sample size (almost 17,000 samples) and the overdispersion observed in earlier analysis of these data without HapMap imputation4. It has been argued that the greater contributor to overdispersion in these data is bias (such as differential genotyping error), rather than population structure4; we therefore carefully examined cluster plots for all SNPs used to impute associated SNPs. Three loci showed suggestive evidence for association (P < 10−7) in regions not previously associated with T1D (Supplementary Fig. 1 and Supplementary Table 2). One SNP, rs229484, was proximal (30 kb) to a nearby known T1D locus (most associated SNP, rs229541)3, also at 22q13.1, but was separated by two moderate recombination hotspots, and there was low linkage disequilibrium (LD) between the two markers (r2 = 0.1; D′ = 0.4).

To replicate these potential effects, we carried out direct genotyping of the three SNPs using TaqMan in a subset of the GWA samples, additional case-control and family samples, and we obtained evidence for association in two of the three loci (Table 1 and Supplementary Table 3). In these two loci, the overall levels of significance were <10−8 (P = 4.13 × 10−9 for rs2304256 and P = 1.62 × 10−10 for rs941576).

Table 1 Association testing of two SNPs using direct genotyping in case-control and family samples
Full size table

rs2304256:C>A (OR for A versus C = 0.86) is located within the TYK2 gene at chromosome 19p13.2, which is implicated in interferon-α, interleukin (IL)-6, IL-10 and IL-12 signaling. This is a region of wide LD containing several functional candidate genes (Supplementary Fig. 2). rs2304256 is one of six SNPs in the 1000 Genomes Project database (pilot 1, April 2009) in mutual tight LD (r2 > 0.9); two are located within TYK2 (rs34725611 and rs11085725 in introns 6 and 23, respectively) and the remaining three (not yet in dbSNP) are downstream of TYK2 and upstream of ICAM3. No other SNPs had r2 > 0.62 with any of these six SNPs. There is evidence that TYK2 is involved in multiple autoimmune diseases: a low-frequency nonsynonymous TYK2 SNP (rs34536443:G>C, P1104A, minor allele frequency 0.04 versus 0.29 for rs2304256) has been convincingly associated with multiple sclerosis6 and, in smaller samples and at lower statistical significance, with ankylosing spondylitis and autoimmune thyroid disease8. The T1D associated rs2304256 is itself a nonsynonymous SNP (V362F) that has also been associated with systemic lupus erythematosus5. In all five diseases, the minor—and inferred nonancestral9—allele (A, encoding phenylalanine, for rs2304256 and T1D and lupus; C, encoding alanine, for rs34536443 and multiple sclerosis, ankylosing spondylitis and autoimmune thyroid disease)5,6,8.

Notably, the newly identified locus with the strongest association with T1D susceptibility occurred in a well-established imprinted region on chromosome 14q32.2 (ref. 10), marked by SNP rs941576:A>G (OR for G versus A = 0.9). Beyond the insulin T1D susceptibility locus, marked by rs7111341 (ref. 4), we do not know of any other T1D SNPs in established imprinted genes. Within this imprinted region of just over 1 Mb, a mixture of paternally derived (DLK1, RTL1 and DIO3) and maternally derived (BEGAIN, MEG3, MEG8 and DIO3OS) genes are expressed10 (Fig. 1). We therefore tested for a parent-of-origin effect, expecting to see excess transmissions of the risk allele from either fathers or mothers (but not both) if the SNP is influencing one of these imprinted genes. A simple way to do this is to consider separately the paternal and maternal transmissions in a transmission disequilibrium testing framework; this revealed strong evidence for reduced paternal transmission of the protective G allele (P = 6.3 × 10−8). Although the maternal transmissions are distorted in the same direction, and a small effect of the maternal copy cannot be discounted, there was no significant evidence for such an effect (P = 0.11; Table 2).

Figure 1: Imprinted region on chromosome 14q32.2.
figure 1

Region shown is delimited by most distant genes known to be imprinted10 with positions according to the NCBI36 assembly of the human genome. Top panel shows –log10(P) from 1–degree of freedom tests of association with SNPs across the region. Black, SNPs directly genotyped; blue, SNPs imputed from HapMap. Middle panel shows location and orientation of genes in the region. Blue, paternally expressed genes; black, maternally expressed genes. Bottom panel shows recombination rates (cM/Mb) from HapMap. Solid green line indicates location of rs941576 in all panels.

Full size image
Table 2 Transmission disequilibrium tests of rs941576:A>G
Full size table

Effects resulting from the action of maternal genotype in utero are confounded with imprinting effects11, so we fitted a model allowing for both maternal genotype and imprinting effects. This has been approached in case-parent trio data by log-linear modeling of counts of trios by parental and affected offspring genotype. We extended this method to allow for the fact that many of our families had multiple affected offspring (see Online Methods). The imprinting-only model was preferred (Supplementary Table 4); under that model, the imprinting effect was highly significant (P = 1.85 × 10−8; ratio of allelic effects for paternally to maternally inherited alleles = 0.75). This test gains power by using information on parental asymmetry induced by parent-of-origin effects. Asymmetry was clearly shown in our data: the protective allele (G) was less common among fathers of affected offspring than among mothers (0.43 versus 0.47, respectively; P = 6.53 × 10−7). To confirm that the results were not falsely positive, driven by unusual patterns in a subset of the data, we reanalyzed the families subdivided by broad geographical region and found consistent effect estimates across all regions (Table 3).

Table 3 Imprinting analysis of rs941576:A>G
Full size table

The SNP rs941576 lies within intron 6 of the maternally expressed noncoding RNA gene MEG3. However, our observation that only paternal transmissions alter T1D risk suggests that the causal variant influences one of the paternally expressed imprinted genes in its neighborhood (DLK1, RTL1 or DIO3). rs941576 is between and downstream of both DLK1 and RTL1 and upstream of DIO3, at distances of 105 kb, 41 kb and 721 kb respectively. Unusually for a locus identified from GWA data, the signal is restricted to rs941576, and there are no SNPs in HapMap or the current prerelease of the 1000 Genomes Project (pilot 1, April 2009) that are in strong or moderate LD with rs941576 (all r2 < 0.5; data not shown). Although that does not preclude the existence of an as-yet-unknown variant (SNP or structural variant) in tighter LD, rs941576 lies within a region conserved across mammalian species, including opossum. This is notable because the region in opossum is not imprinted and shows no sequence homology to MEG3, and although it shows some sequence homology to mouse Rtl1 and human RTL1, the opossum Rtl1 sequence seems to be extensively degraded12. Thus, if the region is conserved because it contains regulatory elements of nearby genes, these must regulate one of the genes common to all mammals (DLK1 or DIO3).

Although rs941576 lies some distance from the paternally expressed genes in the region, regulatory regions can lie >100 kb from their target genes, particularly in imprinted regions13. This region is already subject to long-range cis-acting regulation from the intergenic differentially methylated region located 12.5 kb upstream of MEG3 (ref. 14). Insertion of a transgene in the mouse downstream of this differentially methylated region causes loss of imprinting on the paternal chromosome, biallelic expression of Meg3 (previously known as Gtl2) and reduced expression of Dlk1 (ref. 15). Thus, it is plausible that this SNP (or another unknown variant nearby) alters the regulation of the paternally expressed DLK1 or RTL1.

Of the paternally expressed genes, only DLK1 has a strong functional candidacy. It is most strongly expressed in human heart, pancreatic islet cells, pituitary tissue, ovaries, placenta and testes (T1DBase, BioGPS), is related to members of the Notch-Delta family of signaling molecules and encodes a membrane-bound protein that can be cleaved to form fetal antigen-1 (FA1)16. FA1 is involved in differentiation of many cell types17, including pancreatic beta cells, where FA1 immunoreactivity has been localized to glucagon-negative cells in the mature pancreas18. FA1 is also involved in hematopoiesis, including differentiation and function of B lymphocytes19,20, and has been shown to increase expression of proinflammatory cytokines in human bone marrow mesenchymal stem cells and promote B-cell proliferation in human peripheral blood21. Thus, there are a number of ways in which variation in DLK1 expression could alter susceptibility to T1D, which is caused by autoimmune destruction of insulin-producing beta cells in the pancreas.

The mechanisms underlying imprinting are not fully understood but are known to involve epigenetic processes, including DNA methylation and histone acetylation. The causal variant underlying this association could directly alter the expression of the paternally inherited copy of a nearby gene (DLK1 seems to be the strongest candidate), or it could interfere subtly with the imprinting mechanism and in turn alter expression of either the paternally or maternally inherited copies of a target gene. Although rs941576 may be tagging an unknown causal variant, there is support for the hypothesis that this SNP is itself the causal variant, given its isolation from other SNPs in terms of LD and its location in a conserved and presumably regulatory region.

Rare disorders related to imprinting defects are known (such as Prader-Willi syndrome, MIM#176270). For common complex diseases, over 300 reproducibly associated1 loci have been reported, but we are not aware of any convincing evidence for another susceptibility locus subject to parent-of-origin effects. At least one common disease locus overlaps a known imprinted region: the T1D-associated region of chromosome 11p15 contains the genes encoding insulin and IGF2, but a previous report by our group of potential parent-of-origin effects at this locus in T1D22 has not yet been substantiated. We are aware of only one other report of a parent-of-origin effect, in basal cell carcinoma23, although this was only shown in a single population and at a relatively modest level of statistical significance (P = ∼0.01).

Methods

Sample selection and genotyping.

A total of 7,514 cases and 9,045 control samples were included from three GWA studies: WTCCC (UK), T1DGC (UK) and GoKinD-NIMH (USA). The samples and their genotyping have been described2,3,4. Numbers of samples from each study and genotyping platform are given in Supplementary Table 1. SNP and sample-exclusion criteria were as applied previously4. Briefly, all subjects were of self-reported white European ancestry; samples were excluded if they showed evidence of non-European ancestry, or if they duplicated or were closely related to another sample in the study. SNPs were excluded if the minor allele frequency fell below 1% in cases or controls, if they deviated from Hardy-Weinberg equilibrium (P < 5.7 × 10−7), if the call rate fell below 95% (WTCCC and T1DGC) or if a genotype-calling metric indicated insufficient separation of the signal clouds (GoKinD-NIMH)24.

SNPs showing suggestive association in the imputed analysis were genotyped directly using TaqMan (Applied Biosystems) on a subset of the GWA samples (the T1DGC, all WTCCC cases and about half the WTCCC controls were available to us), additional case-control samples and a set of family samples with T1D-affected offspring (Supplementary Table 1). The additional case and control samples have also been described4. The family samples were drawn from across Europe and America and were predominantly of self-reported white European origin; we did not exclude subjects who self-reported a nonwhite European origin, as testing for transmissions within families is equivalent to a pseudo–case-control approach with ethnically matched controls. All TaqMan genotyping data were scored twice to minimize error; the second operator was unaware of case-control status and family structure.

Imputation.

For each of the three GWA studies, we divided SNPs from HapMap version 2 (release 24) into two sets: those that were genotyped and passed quality control thresholds in the study (X), and those that were not genotyped or failed quality control (Y). The R package snpMatrix7 from the BioConductor project25 was to used calculate imputation ′rules′ for prediction of each SNP in Y from nearby SNPs in X using HapMap genotypes and to carry out association tests for the imputed SNPs. The algorithms used in snpMatrix, together with the parameter settings we used, are described below.

In regions of high LD, the genotype of one SNP can be related to the genotypes of others by a linear regression26,27,28. The first step in calculating an imputation rule is to select a set of ′tag′ SNPs by forward stepwise regression of the Y SNP on the nearest 50 X SNPs (subject to a maximum missing-data requirement). New SNPs are added to the regression until either (i) R2 > 0.95, (ii) the change in R2 is <0.05 or (iii) the number of tag SNPs reaches four. Regression calculations are carried out at the genotype level, with each SNP genotype coded 0, 1 or 2. If a prediction of R2 ≥ 0.95 cannot be achieved using this stepwise regression approach, then an alternative imputation rule is attempted using the set of tag SNPs selected by the forward stepwise procedure. Using the conventional expectation maximization algorithm, frequencies are estimated for the haplotypes of the Y SNP plus the selected tags. Conditional probabilities of the Y allele given the tag SNP haplotype are calculated and provide the imputation rule. This rule is used in preference to the regression rule if the improvement in R2 exceeds 0.1.

These imputation rules are then applied to the main study data set to calculate the expectation of each Y SNP conditional on typed SNPs. This expectation is not generally an integer, and the Cochran-Armitage test then becomes a t test comparing the mean imputation score in cases with that in controls. Extension to allow for stratified comparisons and to combine information from different studies is straightforward: differences between mean scores are simply averaged over strata (and studies), with weights inversely proportional to their variances. These procedures are all implemented in snpMatrix.

This imputation method is computationally faster than those based on hidden Markov models29 or variable-length Markov chains30. For a subset of our data, we compared our imputation results with those from IMPUTE29 and found them to be very similar. It has an additional advantage over such methods in that, because each imputation is based on a small number of tag SNPs, it is easier to differentiate between genuine associations and those caused by poor clustering and differential measurement error; for each putative association, allele signal plots for all tags were visually inspected.

Association analysis.

Single SNP association score tests were conducted for each HapMap SNP within each cohort using direct genotypes if available, or imputed genotypes if not. The score is calculated using the equation

where Yi and Xi are the phenotype (case or control) and genotype data, respectively, for subject i. When a SNP is not directly observed, Xi is replaced by its expected value calculated under the null hypothesis as described above. When it is poorly imputed, this expected value is shrunk toward X̄ and contributes little to the test statistic. The permutation variance (the variance under random permutation of Y) is used to calculate the χ2 test. The score statistics were combined first across strata within cohorts and finally across cohorts using the method proposed by Mantel31. The scores (Ui, where i denotes cohort or stratum) and the variances (Vi) are summed to form an overall test of association, (∑ Ui)T (∑ Vi)−1(∑ Ui). Strata were defined by UK region in the case of the WTCCC and T1DGC samples, and by an estimated ancestry score derived from principal components in the case of the GoKinD-NIMH samples3. Testing for association with SNPs on the X chromosome was carried out using a previously proposed method32. Overdispersion of the test statistics was calculated after removal of known T1D loci4, and these parameters were used to calculate the adjusted P values given in Supplementary Table 2.

SNPs showing overall association (P < 1 × 10−7) in regions not previously reported4 were subject to further screening. Cluster plots of each SNP used for imputation were examined manually, and the results were discarded unless all cluster plots for all cohorts were considered clearly separated. One of the cohorts studied (USA) was not designed as a T1D case-control study and was serendipitously assembled after cases and controls were genotyped on different versions of the Affymetrix 500K chip and to different protocols. This cohort was subject to greater differential bias than were the other cohorts. As a result, many SNPs were found that showed (often extreme) association in the USA samples (P < 1 × 10−7) but no association in the T1DGC and WTCCC samples combined (P > 1 × 10−3); for these SNPs, only the data from T1DGC and WTCCC were combined.

Family data were analyzed by transmission disequilibrium testing, splitting multiplex families into parent offspring trios and using a pseudo–case-control framework to estimate allelic effects. A score statistic was also generated, and a score test for association in case-controls and families combined was conducted by summing the scores and variances as described above.

Imprinting test.

We used a logistic regression approach to test for imprinting and maternal genotype effects on risk in offspring. This approach was originally proposed by Weinberg10,33 for data consisting of trios of an affected individual and both parents, but we required extension to deal with our data, which included families with multiple affected offspring. Weinberg's approach is to analyze counts of case-parent trios classified by genotype of mother (M), father (P) and affected offspring (O) in a 3 × 3 × 3 table. Of the 15 cells in this table consistent with Mendelian transmission, five concern families in which the genotypes of the two parents are concordant; these are not informative in the analysis. The remaining ten cells can be organized by mating type and offspring genotype into five pairs in which the maternal and paternal genotypes are considered interchangeable (Supplementary Table 5). In the absence of maternal genotype and imprinting effects, and assuming that, in the population from which families are drawn, the two possible parental genotype combinations within each mating type are equally frequent, their frequencies in case-parent trios will also not differ systematically. However, maternal genotype and imprinting effects will distort these ratios. In Supplementary Table 5, pairs of genotype configurations are set out with the configuration in which the mother carries more copies of the ′2′ allele than the father appearing first. The table also sets out the predictions of a multiplicative model for relative risk conditional on genotype and on parents; the genotype relative risk for the offspring (γ1/1, γ1/2 and γ2/2) is modified by multiplicative effects of the maternal genotype (ϕ1/2 and ϕ2/2, ϕ1/1 being taken as 1) and by a factor θ if a '2' allele was received from the mother rather than from the father. The ratio of these two risks for each mating type gives the ratio of expected frequencies in case-parent trios. This model can be fitted to the observed pairs of case-parent trio frequencies using any standard logistic regression program, thus allowing estimation and testing of maternal genotype and imprinting effects.

Extension of this method to deal with families in which there may be several affected offspring is relatively straightforward. Again, we tabulated counts of families by genotype of mothers, fathers and offspring, but there were then more possible cells in the tabulation. For example, with two affected offspring, there are seven informative pairs of genotype configurations (Supplementary Table 6). Under the assumption that the SNP under observation is the sole causal variant or has r2 = 1 with a sole causal variant, disease occurrences in the offspring are conditionally independent given their genotypes and their parents, and the ratio of expected frequencies is given by the ratio of products of predicted relative risks for the two offspring. Extension to the case of more than two affected offspring follows similar principles. For families with three affected offspring, there are 9 informative pairs of genotype configuration, for four affected offspring, 11, and so on. Logistic regression can then be used to estimate and test for effects of maternal genotype and imprinting in the general case, as in our study, where the data consist of families with varying numbers of affected offspring.

In the case where the SNP tested is not the sole causal variant (or in perfect LD with it), disease occurrences in offspring are not conditionally independent and there may be some bias. We would expect this to be small when the SNP has high r2 with the causal variant. Moreover, the type 1 error rate will be unaffected by departure from conditional independence when testing the hypothesis of no imprinting and no maternal genotype effect against presence of either (or both) effects, although the method may then not be fully efficient.

URLs.

1000 Genomes, http://www.1000genomes.org; BioGPS, http://biogps.gnf.org; International HapMap Project, http://www.hapmap.org; T1DBase, http://www.t1dbase.org.