The imprinted DLK1-MEG3 gene region on chromosome 14q32.2 alters susceptibility to type 1 diabetes.

Genome-wide association (GWA) studies to map common disease susceptibility loci have been hugely successful, with over 300 reproducibly associated loci reported to date. However, these studies have not yet provided convincing evidence for any susceptibility locus subject to parent-of-origin effects. Using imputation to extend existing GWA datasets, we have obtained robust evidence at rs941576 for paternally inherited risk of type 1 diabetes (T1D; ratio of allelic effects for paternal versus maternal transmissions = 0.75; 95% confidence interval (CI) = 0.71-0.79). This marker is in the imprinted region of chromosome 14q32.2, which contains the functional candidate gene DLK1. Our meta-analysis also provided support at genome-wide significance for a T1D locus at chromosome 19p13.2. The highest association was at marker rs2304256 (odds ratio (OR) = 0.86; 95%CI = 0.82-0.90) in the TYK2 gene, which has previously been associated with systemic lupus erythematosus and multiple sclerosis.

We used imputation to assess association with T1D across 2.6 million polymorphic SNPs from the International HapMap Project in a total of 7514 cases and 9405 controls of European ancestry from three existing genomewide association studies: WTCCC (UK)2, GAIN/NIMH (USA)3, T1DGC (UK)4 (supplementary table 1). The R package snpMatrix6 was used to conduct the imputation and calculate single SNP association score tests for each HapMap SNP. The score tests are 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 an estimated ancestry score derived from principal components in the case of the GoKinD/NIMH samples3). For imputed SNPs, the score statistics are calculated using the expected value of the imputed SNP, given observed SNPs, with the expectation calculated under the null hypothesis.

Overall, there was some over-dispersion of test statistics (λ = 1.14 and λ = 1.09 for 1 degree of freedom (df) and 2df tests respectively). This is consistent with the large sample size (almost 17,000 samples) and the over-dispersion observed in earlier analysis of these data without HapMap imputation4. Barrett et al argue that the greater contributor to over-dispersion in these data is bias (eg differential genotyping error) rather than population structure,4 and therefore cluster plots for all SNPs used to impute associated SNPs were examined carefully. Three loci showed suggestive evidence for association (p < 10−7) in regions not previously associated with T1D (supplementary figure 1​;​ supplementary table 2). One SNP, rs229484, is proximal (30kb) to a nearby known T1D locus (rs2295413), also at 22q13.1, but is separated by two moderate recombination hotspots and there is low LD between the two markers (r2=0.1, D′=0.4).

In order 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 obtained evidence for association in two of the three loci (table 1, supplementary table 3). In these two loci, the overall levels of significance were < 10−8: rs2304256 p = 4.13 × 10−9, rs941576 p = 1.62 × 10−10.

Table 1

Association testing of two SNPs using direct genotyping in case-control and family samples

a. rs2304256 C>A on chromosome 19p13.2
Cohort	N	Fq (A)	Odds ratio (A:C)	(95% CI)	p value
WTCCC	1766/1384	0.299	0.84	(0.75–0.94)	2.68 × 10−3
T1DGC	3838/3883	0.294	0.85	(0.80–0.92)	1.45 × 10−5
Additional	2686/4794	0.290	0.87	(0.81–0.94)	6.02 × 10−4
Families	3099	0.266	0.96	(0.90–1.03)	0.290
Case-control combined	8290/10061	0.293	0.86	(0.82–0.90)	1.43 × 10−10
Families & case-control	(see above)	-	-	-	4.13 × 10−9

Association testing using observed (not imputed) genotypes in a subset of GWA samples, additional case control samples and family samples. SNP names are followed by alleles, ordered as major>minor. N is number of cases/controls, or number of informative transmissions. Fq is the frequency of the minor allele in controls or parents.

rs2304256 C>A (OR for A vs C = 0.86) is located within the TYK2 gene at chromosome 19p13.2, which is implicated in IFN-α, IL-6, IL-10 and IL-12 signalling. This is a region of wide LD containing several functional candidate genes (supplementary figure 2). rs2304256 is one of six SNPs in 1000 Genomes (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. rs2304256 itself is a non-synonymous SNP (Val362Phe) which has also been associated with systemic lupus erythematosus (SLE)5; in both T1D and SLE the minor (and inferred non-ancestral7) allele (A/Phe) appears protective5.

Most interestingly, the newly identified locus with the strongest association with T1D susceptibility occured in a well established imprinted region on chromosome 14q32.28 marked by SNP rs941576 A>G (OR for G vs A = 0.9). Beyond the insulin T1D susceptibility locus, marked by rs7111341 in Barrett et al,4 we do not know of any other T1D SNPs in established imprinted genes. Within this imprinted region of just over 1Mb, a mixture of paternally derived (DLK1, RTL1, DIO3) and maternally derived (MEG3, MEG8) genes are expressed8 (figure 1). Therefore, we 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 was acting to influence 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 (TDT) framework, and this showed 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 is no significant evidence for such an effect (p = 0.11; table 2). However, effects due to the action of maternal genotype in utero are confounded with imprinting effects9, 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 modelling 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 supplementary methods) and found that the imprinting-only model was preferred (supplementary table 4); under that model, the imprinting effect was highly significant (p = 1.85 × 10−8) with the ratio of allelic effects for paternally to maternally inherited alleles equal to 0.75. This test gains power by using information on parental asymmetry induced by parent-of-origin effects. Asymmetry was clearly exhibited in our data: the protective allele (G) is less common amongst fathers of affected offspring than mothers (0.43 vs 0.47, p = 6.53 × 10−7). To reassure ourselves against a false positive result, driven by unusual patterns in a subset of the data, we reanalysed the families subdivided by broad geographical region, and found consistent effect estimates across all regions (table 3).

Figure 1

The imprinted region on chromosome 14q32.2. The region shown is delimited by the most distant genes known to be imprinted 8 with positions according to Hs_NCBI36. The top panel shows -log 10(p) from 1 degree of freedom tests of association with SNPs across the region. SNPs which were directly genotyped are in black, SNPs imputed from HapMap in blue. The second panel shows the location and orientation of genes in the region. Paternally expressed genes are shown in blue, maternally expressed genes in black. The third panel shows recombination rates (cm/Mb) from HapMap. A solid green line shows the location of rs941576 in all panels for reference.

Table 2

Transmission Disequilibrium Tests of rs941576 A>G

Transmissions from	Fq	G Untransmitted	G Transmitted	p value
All parents	0.45	2166	1891	1.6 × 10−5
Fathers	0.43	869	657	6.3 × 10−8
Mothers	0.47	793	730	0.11

Parental frequency (Fq) and transmissions of the rs941576 protective G allele, overall and separated by parent of origin. Frequencies are calculated using all parents. Note that because only transmissions from heterozygous (informative) parents are shown, transmission of a G allele implies non-transmission of A (and vice versa). The sum of maternal and paternal transmissions is less than the number of transmissions from all parents because it is not always possible to identify which parent transmitted which allele.

Table 3

Imprinting analysis of rs941576 A>G

Region	N	exp(-θ̂)	95% CI	p
UK	361	0.792	0.724 – 0.866	9.40 × 10−3
Asia-Pacific	32	0.88	0.656 – 1.18	0.662
Other Europe	257	0.725	0.644 – 0.815	6.08 × 10−3
USA	184	0.764	0.676 – 0.863	0.028
Finland	397	0.697	0.632 – 0.769	2.25 × 10−4
Overall	1231	0.749	0.712 – 0.789	1.85 × 10−8

Imprinting analysis using family data divided by broad geographical region. N is the number of informative families (which is less than the total number of families available, as only transmissions from asymmetric parents are informative). exp(-θ̂) is the ratio of the allelic effect for a paternally inherited risk allele compared to a maternally inherited allele.

The SNP rs941576 lies within intron 6 of the maternally expressed non-coding RNA gene, MEG3. However, our observation that only transmissions from fathers alter T1D risk suggest the causal variant influences one of the paternally expressed imprinted genes in its neighbourhood: DLK1, RTL1 or DIO3. rs941576 is between and downstream of both DLK1 and RTL1 and upstream of DIO3, at distances of 105kb, 41kb and 721kb 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 pre-release of 1000 Genomes Project (pilot 1, April 2009) which 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 interesting because the region is not imprinted in the opossum, there is no sequence homology to MEG3 and, while there is some sequence homology to mouse and human RTL1 gene, it appears to be extensively degraded in opossum10. 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, ie DLK1 or DIO3.

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

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 signalling molecules and encodes a membrane bound protein, which can be cleaved to form fetal antigen 1 (FA1).14 FA1 is involved in differentiation of many cell types15 including pancreatic beta cells where FA1 immunoreactivity has been localised to glucagon-negative cells in the mature pancreas.16 FA1 is also involved in hematopoiesis including differentiation and function of B lymphocytes17,18 and has been shown to increase expression of pro-inflammatory cytokines in human bone marrow mesenchymal stem cells and promote B cell proliferation in human peripheral blood.19 Thus there are a number of ways in which variation in the expression of DLK1 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 yet fully understood, but are known to involve epigenetic processes including DNA methylation and histone acetylation. The causal variant underlying this association could be acting directly to alter the expression of the paternally inherited copy of a nearby gene (DLK1 appears to be the strongest candidate), or it could act by interfering 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 linkage disequilibrium, and its location in a conserved and, presumably, regulatory region.

Rare disorders related to imprinting defects are known (eg Prader-Willi syndrome, OMIM 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 insulin and IGF2 genes, but a previous report by our group of potential parent of origin effects at this locus in T1D20 has not yet been substantiated. We are aware of only one other report of a parent of origin effect, in basal cell carcinoma,21 although this was only demonstrated in a single population and at a relatively modest level of statistical significance (p ≈ 0.01).