Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes.

Type 1 diabetes (T1D) is a common autoimmune disorder that arises from the action of multiple genetic and environmental risk factors. We report the findings of a genome-wide association study of T1D, combined in a meta-analysis with two previously published studies. The total sample set included 7,514 cases and 9,045 reference samples. Forty-one distinct genomic locations provided evidence for association with T1D in the meta-analysis (P < 10(-6)). After excluding previously reported associations, we further tested 27 regions in an independent set of 4,267 cases, 4,463 controls and 2,319 affected sib-pair (ASP) families. Of these, 18 regions were replicated (P < 0.01; overall P < 5 × 10(-8)) and 4 additional regions provided nominal evidence of replication (P < 0.05). The many new candidate genes suggested by these results include IL10, IL19, IL20, GLIS3, CD69 and IL27.

Results from linkage and association studies in T1D have long supported a model in which the major risk factor for T1D resided in the HLA region on chromosome 6p21. Candidate gene studies carried out over a number of years identified four non-HLA T1D risk loci: INS, CTLA4, PTPN22, and IL2RA1-4. Recently, the application of genome-wide SNP typing technology to large sample sets and comparisons with results from other immune-mediated diseases have provided convincing support for 19 additional T1D loci5-13, all with allelic odds ratios (OR’s) of less than 1.3.

In order to have adequate power to detect additional T1D risk loci with ORs in the range of 1.1 to 1.3, we performed a new genome-wide association scan using British cases and controls and used this dataset in a meta-analysis which included 7,514 cases and 9,045 reference samples (Table 1). The other datasets included in the meta-analysis were from the Wellcome Trust Case Control Consortium (WTCCC) study7 and a study12 that utilized T1D cases from the Genetics of Kidneys in Diabetes (GoKinD) study of diabetic nephropathy14, and reference samples from the National Institute of Mental Health (NIMH) study15.

Table 1

Samples from three genome-wide association analyses of type 1 diabetes used in this analysis.

Subjectsa	GWA Meta-analysis				Replication Study				Total
	T1DGC	GoKinD/NIMH	WTCCC	Combined	T1DGC	UK	Danish	Combined
Cases	3,983	1,601	1,930	7,514	-	2,499	1,768	4,267	11,781
Reference	3,999	1,704	3,342	9,045	-	2,690	1,980	4,670	13,715
Totals	7,982	3,305	5,272	16,559	-	5,189	3,748	8,937	25,496
Triosb	-	-	-	-	4,342	-	-	-	4,342

The derivation of subjects from the various indicated studies is described in detail in the Methods section.

From 2,319 affected sib-pair families.

The two earlier studies (WTCCC and GoKinD/NIMH) used Affymetrix 500K platforms while the new (T1DGC) study used the Illumina 550K platform. Of the 841,622 SNPs genotyped in these studies which had minor allele frequencies (MAF) exceeding 1% and passed our quality control standards, 328,044 were only genotyped by the Affymetrix platform, 437,739 only by the Illumina platform, and 75,839 were genotyped by both platforms. Since only 9% of SNPs are shared between these platforms, imputation was used to combine results across studies. To develop imputation rules, we took advantage of the fact that 1,422 of the original WTCCC controls which were included in the T1DGC study had been genotyped on both platforms (Methods).

An analysis using Mantel’s extension to the 1 degree of freedom (1 df) Cochran-Armitage trend test which combined comparisons over the three studies yielded 41 distinct genomic locations with P-values < 10-6 (Figure 1) (Individual plots for each study are in Supplementary Figure 1). Fifteen of these sites were in regions where there were prior reports of association to T1D (Table 2). The remaining 26 of these locations along with one weaker association on the X chromosome, were chosen for further analysis. To address the possible effects of population structure, the analyses were stratified by geographical region in the case of the British studies and by a “propensity score” based on principal components analysis on the US study. This was only partially successful in reducing the over-dispersion of test statistics, a large part of which derived from the US data (Table 3). If the residual over-dispersion were due to population structure, there would be a strong case for correcting the P-values (as shown in Table 3). However, the modest effect of the stratified analysis on over-dispersion, taken together with the absence of any over-dispersion in case-only interaction tests (see below) suggests that it is more likely due to differential genotyping errors. In this case, correction of the most significant P-values would be over-conservative since we have carefully checked all genotyping cluster plots for associated SNPs. The genomic control corrected P-values are nevertheless shown in Supplementary Table 1. The strongest associations tended to become somewhat less significant, but the choice of regions for follow-up, based on the criteria of P < 10-6, was not affected. We also carried out, for SNPs with minor allele frequency exceeding 10%, 2 df “genotype” tests which would be more sensitive to associations showing marked dominance (deviation from an additive model, on the log scale). Significance was notably increased, by 3 to 4 orders of magnitude, at three SNPs, but was less significant than the corresponding 1 df tests otherwise (Supplementary Table 1) yielding no additional findings at P < 10-6. The results of both simple and stratified 1 df tests of these SNPs, separated by study, are shown in Supplementary Tables 3 and 4. Quantile-quantile plots for tests in our new (T1DGC) study, and in the meta-analysis, after removal of tests for SNPs in linkage disequilibrium (LD) regions surrounding known and putative associations, are shown in Supplementary Figure 2a and 2b.

Figure 1

Genome-wide plots of -log10 P-values from stratified 1 df tests combining results from all three studies. Values of -log10 P greater than 10 are plotted at 10. SNPs only present on the Illumina chip are shown in blue, those only present on the Affymetrix chip in red, and those present on both chips are shown in black. Points are plotted in the order red, blue, black. Previously known disease susceptibility loci are marked by vertical black lines, while new findings from the current analysis are marked by vertical grey lines (solid lines for convincingly replicated loci and dashed lines for nominally replicated results).

Table 2

Results for locations of known susceptibility loci for type 1 diabetes.

SNPa	Chromosome	LD region	GWA p-value	Gene of Interestc	References
rs2476601	1p13.2	113.62-114.46	8.5 × 10-85	PTPN22	2
rs2816316	1q31.2	190.73-190.82	3.1 × 10-5	RGS1	10
rs917997	2q12.1	102.22-102.58	0.067b	IL18RAP	10
rs1990760	2q24.2	162.67-163.10	6.6 × 10-9	IFIH1	5
rs3087243	2q33.2	204.38-204.53	1.2 × 10-15	CTLA4	3
rs11711054	3p21.31	45.96-46.63	1.7 × 10-5	CCR5	10
rs4505848	4q27	123.13-123.83	4.7 × 10-13	IL2	9,12
rs6897932	5p13.2	35.84-36.07	0.026	IL7R	9
rs9268645	6p21.32	24.70-34.00	<< 10-100	MHC	16
rs11755527	6q15	90.86-91.10	5.4 × 10-8	BACH2	12
rs2327832	6q23.3	137.80-138.40	0.0003	TNFAIP3	11
rs1738074	6q25.3	159.13-159.62	0.006	TAGAP	10
rs12251307	10p15.1	6.07-6.24	1.3 × 10-13	IL2RA	4
rs11258747	10p15.1	6.48-6.59	1.2 × 10-7	PRKCQ	12
rs7111341	11p15.5	2.02-2.26	4.4 × 10-48	INS	1
rs2292239	12q13.2	54.64-55.09	2.2×10-25	ERBB3	9, 30
rs3809114	12q13.3	55.23-57.27	0.002	multiple	13
rs3184504	12q24.12	109.77-111.72	2.8 × 10-27	SH2B3	9
rs3825932	15q25.1	76.77-77.05	7.7 × 10-8	CTSH	12
rs12708716	16p13.13	10.92-11.56	2.2×10-16	CLEC16A	8, 9
rs1893217	18p11.21	12.73-12.92	3.6 × 10-15	PTPN2	9
rs763361	18q22.2	65.63-65.72	1.2 × 10-5	CD226	9
rs11203203	21q22.3	42.68-42.76	1.7×10-9	UBASH3A	6
rs229541	22q13.1	35.90-36.00	2.1 × 10-7	C1QTNF6	12

Focal SNP in each region was taken from the referenced studies.

2d.f. test, as this effect does not conform to a multiplicative model.

The gene of interest choice was based on known expression or function in the immune system, association results from other immune-mediated diseases, the extent of the region of LD based on recombination frequencies from HapMap data, and the location of the SNPs with the highest T1D association; this selection does not infer that this is the causal gene in the region. Other genes, recombination frequency and summary association results are shown in T1DBase.

Table 3

Over-dispersion factors (λ) of 1 df association tests

Study	Simple tests		Stratified tests
	λ	λ	p=10-6	p=10-8
WTCCC	1.077	1.062	2.1 × 10-6	2.7 × 10-8
GoKinD/NIMH	1.196	1.150	5.1 × 10-6	9.1 × 10-8
T1DGC	1.066	1.055	1.9 × 10-6	2.4 × 10-8
GB studiesa	1.105	1.092	3.2 × 10-6	5.0 × 10-8
Combinedb	1.136	1.119	3.8 × 10-6	6.0 × 10-8

For the stratified test λ values, the effect of genomic control correction of p-values of 10-6 and 10-8 are also shown.

Values are shown for each study separately and for meta-analyses of both GB studies (WTCCC and T1DGC)

Values are shown for each study separately and for meta-analyses of all three studies.

The most significantly T1D associated SNPs from each of the 27 novel regions selected for replication were genotyped in a further 4,267 cases and 4,670 controls and in 4,342 trios from 2,319 T1DGC families with multiple affected offspring. Genotype data passed design and quality control criteria for 25 of these SNPs. Eighteen regions replicated with P < 0.01 and showed genome-wide significant (P < 5 × 10-8) association in the joint analysis of the genome scans and replication samples (Table 4, individual scan data in Supplementary Table 2). A further three of the remaining seven SNPs also showed P < 0.01 in the replication studies, and a fourth had P < 0.05, but these failed to reach overall P < 5 × 10-8 (Table 4). This study, therefore, adds 18 T1D risk loci to the existing 24, and provides suggestive support for four more. As expected, nearly all of these loci have OR < 1.2, as larger effects would likely have been discovered in earlier studies. Two of the new associations (10q23 and 16q23) contradict this trend and highlight the disparity between genomic coverage of the older Affymetrix 500K chip and the newer Illumina 550K: these loci do not have a good proxy on the Affymetrix chip, explaining why they were not previously identified despite relatively large effect sizes (OR ∼ 1.3).

Table 4

Replication study of new type 1 diabetes risk loci

SNPa	Chr	LD regionb(Mb)	Gene ofinterest (#)c	P-values			RiskAllele	MAFe	OR (95% CI)f
				GWAd	Replication	Combined			Case-control	Families
rs3024505	1q32.1	204.87-205.12	IL10 (5)	2.2×10-6	0.00015	1.9×10-9	C	0.169	0.84 (0.77-0.91)	0.96 (0.88-1.04)
rs10517086	4p15.2	25.64-25.75	(0)	2.8×10-7	0.00021	4.6×10-10	A	0.299	1.09 (1.02-1.17)	1.09 (1.02-1.16)
rs9388489	6q22.32	126.48-127.46	C6orf173 (1)	5.1×10-8	1.4×10-6	4.2×10-13	G	0.452	1.17 (1.10-1.24)	1.05 (0.99-1.12)
rs7804356	7p15.2	26.62-27.17	(10)	3.3×10-8	0.0051	5.3×10-9	T	0.238	0.88 (0.82-0.94)	0.99 (0.92-1.06)
rs4948088	7p12.1	50.87-51.64	COBL (1)	2.7×10-6	0.0019	4.4×10-8	C	0.047	0.77 (0.67-0.90)	0.93 (0.79-1.10)
rs7020673	9p24.2	4.22-4.31	GLIS3 (1)	1.9×10-9	0.00013	5.4×10-12	G	0.502	0.88 (0.83-0.93)	0.97 (0.91-1.03)
rs10509540	10q23.31	90.00-90.27	C10orf59 (1)	6.9×10-9	4.9×10-24	1.3×10-28	T	0.285	0.75 (0.70-0.80)	0.81 (0.76-0.87)
rs4763879	12p13.31	9.51-9.87	CD69 (6)	2.8×10-7	1.1×10-5	1.9×10-11	A	0.368	1.09 (1.02-1.16)	1.12 (1.05-1.19)
rs1465788	14q24.1	68.24-68.39	(2)	1.4×10-8	1.5×10-5	1.8×10-12	G	0.287	0.86 (0.80-0.91)	0.95 (0.89-1.02)
rs4900384	14q32.2	97.43-97.60	(0)	1.1×10-6	0.00042	3.7×10-9	G	0.288	1.09 (1.02-1.16)	1.08 (1.01-1.16)
rs4788084	16p11.2	28.19-28.94	IL27 (24)	5.2×10-8	8.4×10-7	2.6×10-13	G	0.424	0.86 (0.81-0.91)	0.94 (0.88-1.00)
rs7202877	16q23.1	73.76-74.09	(7)	5.7×10-11	1.2×10-6	3.1×10-15	G	0.096	1.28 (1.17-1.41)	1.09 (0.99-1.20)
rs2290400	17q12	34.63-35.51	ORMDL3 (23)	1.3×10-7	8.2×10-7	5.5×10-13	G	0.495	0.87 (0.82-0.93)	0.92 (0.87-0.98)
rs7221109	17q21.2	35.95-36.13	(3)	9.9×10-10	0.0083	1.3×10-9	C	0.353	0.95 (0.89-1.01)	0.94 (0.88-1.00)
rs425105	19q13.32	51.84-52.02	(5)	1.5×10-7	2.6×10-5	2.7×10-11	A	0.162	0.86 (0.79-0.93)	0.90 (0.82-0.98)
rs2281808	20p13	1.44-1.71	(3)	5.0×10-7	4.8×10-6	1.2×10-11	C	0.362	0.90 (0.84-0.95)	0.90 (0.85-0.96)
rs5753037	22q12.2	28.14-29.00	(14)	1.8×10-14	5.8× 10-5	2.6×10-16	T	0.391	1.10 (1.04-1.17)	1.08 (1.02-1.15)
rs2664170	Xq28	153.48-154.10	(16)	3.0×10-5	5.8×10-5	7.8×10-9	G	0.316	1.16 (1.07-1.24)	1.06 (0.97-1.16)
rs2269241	1p31.3	63.87-63.94	PGM1 (1)	5.9×10-6	0.0069	4.2×10-7	G	0.192	1.10 (1.02-1.18)	1.05 (0.98-1.14)
rs1534422	2p25.1	12.53-12.60	(0)	6.7×10-6	0.025	2.1×10-6	G	0.460	1.08 (1.02-1.15)	1.01 (0.95-1.08)
rs12444268	16p12.3	20.17-20.28	(2)	2.0×10-6	0.0045	1.7×10-7	A	0.295	1.10 (1.03-1.17)	1.04 (0.97-1.11)
rs16956936	17p13.1	7.56-7.66	(2)	3.2×10-6	0.0097	5.3×10-7	C	0.135	0.92 (0.84-1.00)	0.92 (0.83-1.01)

SNPs providing evidence of association at P < 0.05 with T1D in replication study. SNPs showing evidence of replication at P < 0.01 and P < 5 × 10-8 overall are listed by autosome 1-22 and chromosome X (n = 18), followed by those SNPs attaining evidence of association in the replication study at P < 0.01 (n = 3) or 0.05 (n = 1) but failing to reach P < 5 × 10-8 overall.

To define an LD region for a given focal SNP, we extended the region to the left until either 0.1 cM had been traversed or until reaching another SNP with p < 10-6. In the latter case we then set this new SNP as the left bound and repeated the process. The right hand boundary was defined in the same way. However, the boundaries of the region 7p12.1 (50.87-51.64 Mb), were chosen on recombination frequency (T1DBase) and the fat that this larger interval contained all of the COBL gene.

Gene names are shown for regions with a functionally interesting candidate or for regions with only one gene. The total number of genes in each LD region are shown in parentheses.

P-values for stratified 1 degree of freedom tests combining data from all three GWA scans in a meta-analysis.

Minor allele frequency in British controls.

Odds ratio (95% CI, confidence interval). Odds ratios represent the effect of a single copy of the indicated allele within the multiplicative model for allelic effects. For rs2664170, on the X chromosome, the model fitted assumes that relative risks for males reflect those between homozygous females25.

The families utilized for replication were derived from affected sib-pair linkage studies. One consequence of ascertainment on the basis of at least two affected siblings was a high frequency of high risk HLA genotypes16. It has been reported that relative risks for several non-HLA loci are reduced in subjects carrying high risk HLA genotypes17, 18, reflecting deviation from a multiplicative model for joint effects, and this would lead us to expect reduced effect sizes in multiple-case families. Indeed, the results of the replication study were generally less convincing in the family data than in the case-control data reflecting smaller effect sizes in the families. One potential explanation for these different effect sizes lies in possible statistical interaction among risk loci leading to a less-than-multiplicative accumulation of risk in samples (such as those from multiplex families) with a large number of risk variants. This hypothesis is difficult to test because power to detect interaction terms is much less than that to find equivalent sized main effects and is doubly compounded when specific causal variants (rather than tag SNPs from a GWA scan) are not known.

We tested for deviation from the model of multiplicative effects with HLA, on a genome-wide basis, by first calculating predictive risk scores using SNPs in the MHC region on each platform, and testing for association between this score and every other SNP in the remainder of the genome. These tests are “case-only” tests for statistical interaction reflecting variation of allelic relative risks with the level of HLA-attributable risk. As noted earlier, these test statistics did not show the over-dispersion which would have been indicative of population stratification (Supplementary Figure 2c). However, the subset of these tests concerning established T1D susceptibility loci tended to have larger chi-squared values than expected by chance (Supplementary Figure 2d). In the majority of cases (31/45), the interaction tests took the opposite sign from the main effect test, consistent with high MHC risk leading to lower risk for other loci. Of the five interactions which reached P < 0.05, four were of this type (loci near 2q24.2/IFIH1, 1p13.2/PTPN22, 17p13.1 and 2q33.2/CTLA4). We carried out a further test by calculating a T1D risk score using all associated loci excluding the MHC region and testing, in cases only, for correlation between this score and the MHC risk score. We found a weak, but significant (P=0.0007) negative correlation, again indicating that risk from HLA and non-HLA sources accumulates at a rate less than expected based on the model of multiplicative effects, so that there is a general tendency for relative risks for non-HLA loci to be reduced when HLA-related risk is high.

Several of the 18 regions identified here contain genes of possible functional relevance to T1D. These include the region 1q32.1 containing the potent immunoregulatory cytokine genes, IL10, IL19 and IL20. The region of strong LD at 9p24.2 contains only a single gene, GLIS3. Mutations in GLIS3 have been reported in children from three different consanguineous families with permanent neonatal diabetes associated with congenital hypothyroidism and other clinical complications19. The region on 12p13.31 harbors a number of immunoregulatory genes including CD69, which is induced by activation of T cells and functions in thymic egress20. Several other members of the calcium-dependent (C-type) lectin (CLEC) domain family with immune functions also map to this region. Overall, our results provide a rich new source of candidate genes, but until further genotyping, re-sequencing and functional studies are performed, it is not possible to be more specific in regard to which genes might be causal.

Methods

Subjects

The WTCCC study has been described elsewhere7. Cases were recruited from pediatric and adult diabetes clinics at 150 National Health Service Hospitals across Great Britain as part of the Genetic Resource for Investigating Diabetes (GRID) collection (www.childhood-diabetes.org.uk/grid.shtml) of the JDRF/WT DIL9. Half of the controls were drawn from the British 1958 Birth Cohort21 and half from a group of blood donors recruited by the WTCCC in collaboration with the UK Blood Services7. The former group was subsequently genotyped on the Illumina 550K platform and was used as controls in the new T1DGC study reported here. Since the removal of this group from the WTCCC study left it somewhat short of controls, we used a group of 1,868 patients with bipolar disorder as additional reference samples — a group conspicuous in the WTCCC studies in its lack of significant differences from control allele frequencies7.

Our new study added approximately 2,500 new controls from the British 1958 Birth Cohort to the 1,500 described above, and compared these with a new group of approximately 4,000 British cases from the JDRF/WT DIL collection. All cases and controls were resident in Great Britain. To minimize the effects of population structure, the case-control comparisons in the WTCCC and T1DGC studies have been stratified by the 12 regions of Great Britain5,7. Sample exclusions in the genome-wide studies are discussed in Supplementary Methods.

Replication studies were carried out in two groups of cases and control as well as 2,319 affected sib-pair families previously recruited and characterized by the T1DGC6. The British cases were from the JDRF/WT DIL, and the controls were drawn from the British 1958 Birth Cohort, and the UK Blood Service controls of the WTCCC. The second set of cases and controls from Denmark were recruited from a nationwide registry. All cases (49% females) were diagnosed before age 18 years and the mean age at onset 9.02 years. Control subjects were randomly selected from the Inter99 study22.

Genotyping

For the T1DGC study, the 4,000 T1D case and 2,500 control DNA samples were selected based on no prior use in a prior genome wide association study and migration as a high molecular weight band of genomic DNA, ∼23 kb, by electrophoresis on a 0.75% agarose gel. All DNA samples were extracted using a chloroform-based method and quantified in triplicate using Picogreen®. Once selected, the case and control DNA were randomized by columns into a 96 well plate format.

For the T1DGC study, genotyping was performed on the Illumina 550K Infinium platform and, for comparability, all genotypes were re-scored using the ILLUMINUS algorithm23. The WTCCC study used the Affymetrix GeneChip Human Mapping 500K Array set, while the GoKinD/NIMH study used genotype data generated with the Affymetrix Genome-wide Human SNP Array 5.0. The 5.0 array incorporates all of the SNPs on the earlier 500K array but on a single chip along with an additional 420K non-polymorphic probes. Details of the scoring of genotypes may be found in the original publications7, 12. The criteria for discarding some SNPs from the analysis are discussed in Supplementary Methods.

For the replication studies, genotyping was performed in a fully blinded fashion using Taqman assays as previously described9.

Statistical methods

One degree of freedom tests are Cochran-Armitage tests for trend alternatives, extended to pool information across multiple studies or across multiple strata within a single study by the method described by Mantel24. The two degree of freedom tests follow similar principles. Testing for association with SNPs on the X chromosome was carried out using the method proposed by Clayton25. More details are given in Supplementary Methods.

The meta-analysis involved studies that used different platforms, necessitating the use of imputation. Since we had a substantial sample typed on both platforms, we used a simple linear regression approach to imputation26. Details of this, and other methods used in the meta-analysis, are given in Supplementary Methods. Supplementary Figure 3 shows the distribution of the quality of imputation, as measured by the coefficient of determination, R2.

Analysis of the replication case-control studies was carried out in a similar manner, by 1 df comparisons of allele frequencies with Danish and UK studies treated as separate strata. The family study was analyzed by the transmission/disequilibrium test (TDT).

The MHC risk score was derived by an adaption of the lasso approach27 to logistic regression of case/control status versus all SNPs in the MHC region (defined as spanning from 24.7 Mb to 34.0 Mb on chromosome 6). This was applied to the combined Affymetrix data, with a dummy variable in the regression to differentiate WTCCC and GoKinD/NIMH studies and, separately to the T1DGC Illumina data. The coefficients for the selected regression equations are shown in Supplementary Table 3. The degree of risk prediction, as demonstrated by the receiver operating curves (Supplementary Figure 4) was very similar in the three study groups.

A case-only test for statistical interaction between each SNP and MHC risk score was carried out by a 1 df test based on the covariance between MHC risk score and the SNP genotype coded 0, 1 or 2. These tests were stratified within study by geographical region or by principal component score, and information pooled across strata and studies as described above. A 2 df test for association, possibly modified by MHC, was calculated by adding the chi-squared interaction test on 1 df to the 1 df chi-squared statistic for the stratified association test.

The lasso analysis of the MHC risk prediction was carried out in the lasso2 package in the R statistical system28. All the remaining analysis was carried out in the snpMatrix package from the bioConductor project 29.