An Expanded Genome-Wide Association Study of Type 2 Diabetes in Europeans.

To characterize type 2 diabetes (T2D)-associated variation across the allele frequency spectrum, we conducted a meta-analysis of genome-wide association data from 26,676 T2D case and 132,532 control subjects of European ancestry after imputation using the 1000 Genomes multiethnic reference panel. Promising association signals were followed up in additional data sets (of 14,545 or 7,397 T2D case and 38,994 or 71,604 control subjects). We identified 13 novel T2D-associated loci (P < 5 × 10-8), including variants near the GLP2R, GIP, and HLA-DQA1 genes. Our analysis brought the total number of independent T2D associations to 128 distinct signals at 113 loci. Despite substantially increased sample size and more complete coverage of low-frequency variation, all novel associations were driven by common single nucleotide variants. Credible sets of potentially causal variants were generally larger than those based on imputation with earlier reference panels, consistent with resolution of causal signals to common risk haplotypes. Stratification of T2D-associated loci based on T2D-related quantitative trait associations revealed tissue-specific enrichment of regulatory annotations in pancreatic islet enhancers for loci influencing insulin secretion and in adipocytes, monocytes, and hepatocytes for insulin action-associated loci. These findings highlight the predominant role played by common variants of modest effect and the diversity of biological mechanisms influencing T2D pathophysiology.

Introduction

Type 2 diabetes (T2D) has rapidly increased in prevalence in recent years and represents a major component of the global disease burden (1). Previous efforts to use genome-wide association studies (GWAS) to characterize the genetic component of T2D risk have largely focused on common variants (minor allele frequency [MAF] >5%). These studies have identified close to 100 loci, almost all of them currently defined by common alleles associated with modest (typically 5–20%) increases in T2D risk (2–6). Direct sequencing of whole genomes or exomes offers the most comprehensive approach for extending discovery efforts to the detection of low-frequency (0.5% < MAF < 5%) and rare (MAF <0.5%) risk and protective alleles, some of which might have greater impact on individual predisposition. However, extensive sequencing has thus far been limited to relatively small sample sizes (at most, a few thousand cases), restricting power to detect rarer risk alleles even if they are of large effect (7–9). Although evidence of rare variant associations has been detected in some candidate gene studies (10,11), the largest study to date, involving exome sequencing in ∼13,000 subjects, found little trace of rare variant association effects (9).

Here, we implement a complementary strategy that makes use of imputation into existing GWAS samples from the DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium with sequence-based reference panels (12). This strategy allows the detection of common and low-frequency (but not rare) variant associations in extremely large samples (13) and facilitates the fine-mapping of causal variants. We performed a European ancestry meta-analysis of GWAS with 26,676 T2D case and 132,532 control subjects, and we followed up our findings in additional independent European ancestry studies of 14,545 T2D case and 38,994 control subjects genotyped using the Metabochip (4). All contributing studies were imputed against the March 2012 multiethnic 1000 Genomes Project (1000G) reference panel of 1,092 whole-genome–sequenced individuals (12). Our study provides near-complete evaluation of common variants with much improved coverage of low-frequency variants, and the combined sample size considerably exceeds that of the largest previous T2D GWAS meta-analyses in individuals of European ancestry (4). In addition to genetic discovery, we fine-mapped novel and established T2D-associated loci to identify regulatory motifs and cell types enriched for potential causal variants, as well as pathways through which T2D-associated loci increase disease susceptibility.

Research Design and Methods

Research Participants

The DIAGRAM stage 1 meta-analyses comprises 26,676 T2D case and 132,532 control subjects (effective sample size Neff = 72,143 individuals, defined as 4/[(1/Ncases) + (1/Ncontrols)]) from 18 studies genotyped using commercial genome-wide single nucleotide variant (SNV) arrays (Supplementary Table 1). The Metabochip stage 2 follow-up comprises 14,545 T2D case and 38,994 control subjects (Neff = 38,645) from 16 nonoverlapping stage 1 studies (4,14). We performed additional follow-up in 2,796 T2D case and 4,601 control subjects from the European Prospective Investigation into Cancer and Nutrition-InterAct (EPIC-InterAct) study (15) and in 9,747 T2D case and 61,857 control subjects from the Resource for Genetic Epidemiology on Adult Health and Aging (GERA) study (16) (Supplementary Material).

Statistical Analyses

We imputed autosomal and X chromosome SNVs using the all-ancestries 1000G reference panel (1,092 individuals from Africa, Asia, Europe, and the Americas [March 2012 release]) using minimac (17) or IMPUTE2 (18). After imputation, from each study we removed monomorphic variants or those with imputation quality r2-hat < 0.3 (minimac) or proper-info <0.4 (IMPUTE2, SNPTEST). Each study performed T2D association analysis using logistic regression, adjusting for age, sex, and principal components for ancestry, under an additive genetic model. We performed inverse-variance weighted fixed-effect meta-analyses of the 18 stage 1 GWAS (Supplementary Table 1). Fifteen of the 18 studies repeated analyses also adjusting for BMI. SNVs reaching suggestive significance P < 10−5 in the stage 1 meta-analysis were followed up. Novel loci were selected using the threshold for genome-wide significance (P < 5 × 10−8) in the combined stage 1 and stage 2 meta-analysis. For the 23 variants with no proxy (r2 ≥ 0.6) available in Metabochip with 1000G imputation in the fine-mapping regions, the stage 1 result was followed up in EPIC-InterAct and GERA (Neff = 40,637), both imputed to 1000G variant density (Supplementary Material). Summary-level statistics from the stage 1 GWAS meta-analysis are available online at http://diagram-consortium.org/downloads.html.

Approximate Conditional Analysis With GCTA

We performed approximate conditional analysis in the stage 1 sample using GCTA v1.24 (19,20). We analyzed SNVs in the 1-Mb window around each lead variant, conditioning on the lead SNV at each locus (Supplementary Material) (21). We considered loci to contain multiple distinct signals if multiple SNVs reached locus-wide significance (P < 10−5), accounting for the approximate number of variants in each 1-Mb window (14).

Fine-Mapping Analyses Using Credible Set Mapping

To identify 99% credible sets of causal variants for each distinct association signal, we performed fine-mapping for loci at which the lead independent SNV reached P < 5 × 10−4 in the stage 1 meta-analysis. We performed credible set mapping using the T2D stage 1 meta-analysis results to obtain the minimal set of SNVs with cumulative posterior probability >0.99 (Supplementary Material).

Type 1 Diabetes/T2D Discrimination Analysis

Given the overlap between loci previously associated with type 1 diabetes (T1D) and the associated T2D loci, we used an inverse-variance weighted Mendelian randomization approach (22) to test whether this was likely to reflect misclassification of T1D case subjects as individuals with T2D in the current study (Supplementary Material).

Expression Quantitative Trait Locus Analysis

To look for potential biological overlap of T2D lead variants and expression quantitative trait locus (eQTL) variants, we extracted the lead (most significantly associated) eQTL for each tested gene from existing data sets for a range of tissues (Supplementary Material). We concluded that a lead T2D SNV showed evidence of association with gene expression if it was in high linkage disequilibrium (LD) (r2 > 0.8) with the lead eQTL SNV (P < 5 × 10−6).

Hierarchical Clustering of T2D-Related Metabolic Phenotypes

Starting with the T2D-associated SNVs, we obtained T2D-related quantitative trait z scores from published HapMap-based GWAS meta-analysis for the following: fasting glucose, fasting insulin adjusted for BMI, HOMA for β-cell function, and HOMA for insulin resistance (23); 2-h glucose adjusted for BMI (24); proinsulin (25); corrected insulin response (CIR) (26); BMI (27); and HDL cholesterol, LDL cholesterol, total cholesterol, and triglycerides (28). When an association result for an SNV was not available, we used the results for the variant in highest LD and only for variants with r2 > 0.6. We performed clustering of phenotypic effects using z scores for association with T2D risk alleles and standard methods (Supplementary Material) (29).

Functional Annotation and Enrichment Analysis

We tested for enrichment of genomic and epigenomic annotations using chromatin states for 93 cell types (after excluding cancer cell lines) from the National Institutes of Health (NIH) Roadmap Epigenomics Project, as well as binding sites for 165 transcription factors from the Encyclopedia of DNA Elements (ENCODE) project (30) and Pasquali et al. (31). Using fractional logistic regression, we then tested for the effect of variants with each cell type and transcription factor annotation on the variant posterior probabilities (πc) using all variants within 1 Mb of the lead SNV for each distinct association signal from the fine-mapping analyses (Supplementary Material). In each analysis, we considered an annotation significant if it reached a Bonferroni-corrected P < 1.9 × 10−4 (i.e., 0.05/258 annotations).

Pathway Analyses With DEPICT

We used the Data-driven Expression Prioritized Integration for Complex Traits (DEPICT) tool (32) to 1) prioritize genes that may represent promising candidates for T2D pathophysiology and 2) identify reconstituted gene sets that are enriched in genes from associated regions and might be related to T2D biological pathways. As input, we used independent SNVs from the stage 1 meta-analysis SNVs with P < 10−5 and lead variants at established loci (Supplementary Material). For the calculation of empirical enrichment P values, we used 200 sets of SNVs randomly drawn from entire genome within regions matching by gene density; we performed 20 replications for false discovery rate (FDR) estimation. Supplementary tables, supplementary material, and DEPICT analyses are available online at http://diagram-consortium.org/2017_Scott_DIAGRAM_1000G/.

Results

Novel Loci Detected in T2D GWAS and Metabochip-Based Follow-up

The stage 1 GWAS meta-analysis included 26,676 T2D case and 132,532 control subjects and evaluated 12.1 million SNVs, of which 11.8 million were autosomal and 260,000 mapped to the X chromosome. Of these, 3.9 million variants had MAF between 0.5 and 5%, a near fifteen-fold increase in the number of low-frequency variants tested for association compared with previous array-based T2D GWAS meta-analyses (2,4) (Supplementary Table 2). Of the 52 signals showing promising evidence of association (P < 10−5) in stage 1, 29 could be followed up in the stage 2 Metabochip data. In combined stage 1 and stage 2 data, 13 novel loci were detected at genome-wide significance (Table 1, Fig. 1, Supplementary Fig. 1A–, and Supplementary Table 3).

Table 1

Novel loci associated with T2D from the combination of 1000G-imputed GWAS meta-analysis (stage 1) and Metabochip follow-up (stage 2)

Locus name*	Stage 1						Stage 2							Stage 1 + stage 2
	Chr:position	SNV†	EA/ NEA	EAF	OR (95% CI)	P value	Chr:position	SNV‡	r2 with lead SNV	EA/ NEA	EAF	OR (95% CI)	P value	OR (95% CI)¢	P value
ACSL1	4:185708807	rs60780116	T/C	0.84	1.09 (1.06–1.13)	7.38 × 10−8	4:185714289	rs1996546	0.62	G/T	0.86	1.08 (1.03–1.13)	5.60 × 10−4	1.09 (1.06–1.12)	1.98 × 10−10
HLA-DQA1	6:32594309	rs9271774	C/A	0.74	1.10 (1.06–1.14)	3.30 × 10−7	6:32594328	rs9271775	0.91	T/C	0.80	1.08 (1.03–1.13)	7.59 × 10−4	1.09 (1.06–1.12)	1.11 × 10−9
SLC35D3	6:137287702	rs6918311	A/G	0.53	1.07 (1.04–1.10)	6.67 × 10−7	6:137299152	rs4407733	0.92	A/G	0.52	1.05 (1.02–1.08)	1.63 × 10−3	1.06 (1.04–1.08)	6.78 × 10−9
MNX1	7:157027753	rs1182436	C/T	0.80	1.08 (1.05–1.12)	8.30 × 10−7	7:157031407	rs1182397	0.92	G/T	0.85	1.06 (1.02–1.11)	4.38 × 10−3	1.08 (1.05–1.10)	1.71 × 10−8
ABO	9:136155000	rs635634	T/C	0.18	1.08 (1.05–1.12)	3.59 × 10−7	9:136154867	rs495828	0.83	T/G	0.20	1.06 (1.01–1.10)	1.23 × 10−2	1.08 (1.05–1.10)	2.30 × 10−8
PLEKHA1	10:124186714	rs2292626	C/T	0.50	1.09 (1.06–1.11)	1.75 × 10−12	10:124167512	rs2421016	0.99	C/T	0.50	1.05 (1.02–1.08)	2.30 × 10−3	1.07 (1.05–1.09)	1.51 × 10−13
HSD17B12	11:43877934	rs1061810	A/C	0.28	1.08 (1.05–1.11)	5.29 × 10−9	11:43876435	rs3736505	0.92	G/A	0.30	1.05 (1.01–1.08)	4.82 × 10−3	1.07 (1.05–1.09)	3.95 × 10−10
MAP3K11	11:65364385	rs111669836	A/T	0.25	1.07 (1.04–1.10)	7.43 × 10−7	11:65365171	rs11227234	1.00	T/G	0.24	1.05 (1.01–1.08)	8.77 × 10−3	1.06 (1.04–1.09)	4.12 × 10−8
NRXN3	14:79945162	rs10146997	G/A	0.21	1.07 (1.04–1.10)	4.59 × 10−6	14:79939993	rs17109256	0.98	A/G	0.21	1.07 (1.03–1.11)	1.27 × 10−4	1.07 (1.05–1.09)	2.27 × 10−9
CMIP	16:81534790	rs2925979	T/C	0.30	1.08 (1.05–1.10)	2.72 × 10−8	16:81534790	rs2925979	1.00	T/C	0.31	1.05 (1.02–1.08)	3.06 × 10−3	1.07 (1.04–1.09)	2.27 × 10−9
ZZEF1	17:4014384	rs7224685	T/G	0.30	1.07 (1.04–1.10)	2.00 × 10−7	17:3985864	rs8068804	0.95	A/G	0.31	1.07 (1.03–1.11)	4.11 × 10−4	1.07 (1.05–1.09)	3.23 × 10−10
GLP2R	17:9780387	rs78761021	G/A	0.34	1.07 (1.05–1.10)	5.49 × 10−8	17:9791375	rs17676067	0.87	C/T	0.31	1.03 (1.00–1.07)	3.54 × 10−2	1.06 (1.04–1.08)	3.04 × 10−8
GIP	17:46967038	rs79349575	A/T	0.51	1.07 (1.04–1.09)	2.61 × 10−7	17:47005193	rs15563	0.78	G/A	0.54	1.04 (1.01–1.07)	2.09 × 10−2	1.06 (1.03–1.08)	4.43 × 10−8

*The nearest gene is listed; this does not imply this is the biologically relevant gene. †Lead SNV types: all map outside transcripts except rs429358 (missense variant) and rs1061810 (3′ untranslated region). ‡Stage 2: proxy SNV (r2 > 0.6 with stage 1 lead SNV) was used when no stage 1 SNV was available. ¢The meta-analysis OR is aligned to the stage 1 SNV risk allele. Chr, chromosome; EA, effect allele; EAF, effect allele frequency; NEA, noneffect allele.

Figure 1

The effect sizes of the established (blue diamonds, N = 69, P < 5 × 10−4) (Supplementary Material), novel (red diamonds, N = 13), and additional distinct (sky blue diamonds, N = 13) (Supplementary Table 7) signals according to their risk allele frequency (Supplementary Table 3). The additional distinct signals are based on approximate conditional analyses. The distinct signal at TP53INP1 led by rs11786613 (Supplementary Table 7) is plotted (sky blue diamond). This signal did not reach locus-wide significance but was selected for follow-up because of its low frequency and absence of LD with previously reported signal at this locus. The power curve shows the estimated effect size for which we had 80% power to detect associations. Established common variants with OR >1.12 are annotated.

Lead SNVs at all 13 novel loci were common. Although detected here using 1000G imputed data, all 13 were well captured by variants in the HapMap CEU (Central EUrope) reference panel (two directly, 10 via proxies with r2 > 0.8, and one via proxy with r2 = 0.62) (Supplementary Material). At all 13, lead variants defined through 1000G and those seen when the SNP density was restricted to HapMap content had broadly similar evidence of association and were of similar frequency (Supplementary Fig. 2 and Supplementary Table 3). Throughout this article, loci are named for the gene nearest to the lead SNV, unless otherwise specified (Table 1 and Supplementary Material).

Adjustment for BMI revealed no additional genome-wide significant associations for T2D and, at most known and novel loci, there were only minimal differences in statistical significance and estimated T2D effect size between BMI-adjusted and unadjusted models. The four signals at which we observed a significant effect of BMI adjustment (Pheterogeneity <4.4 × 10−4; based on 0.05/113 variants currently or previously reported to be associated with T2D at genome-wide significance) were FTO and MC4R (at which the T2D association is known to reflect a primary effect on BMI) and TCF7L2 and SLC30A8 (at which T2D associations were strengthened after BMI-adjustment) (Supplementary Fig. 3 and Supplementary Table 4).

Insights Into Genetic Architecture of T2D

In this meta-analysis, we tested 3.9 million low-frequency variants (r2 ≥ 0.3 or proper-info ≥0.4; minor allele present in ≥3 studies) for T2D association, constituting 96.7% of the low-frequency variants ascertained by the 1000G European panel (March 2012) (Supplementary Table 2). For variants with risk allele frequencies (RAF) of 0.5%, 1%, or 5%, we had 80% power to detect association (P < 5 × 10−8) for allelic odds ratios (ORs) of 1.80, 1.48, and 1.16, respectively, after accounting for imputation quality (Fig. 1 and Supplementary Table 5). Despite the increased coverage and sample size, we identified no novel low-frequency variants at genome-wide significance (Fig. 1).

Since we had only been able to test 29 of the 52 promising stage 1 signals on the Metabochip, we investigated whether this failure to detect low-frequency variant associations with T2D could be a consequence of selective variant inclusion on the Metabochip. Among the remaining 23 variants, none reached genome-wide significance after aggregating with GWAS data available from EPIC-InterAct. Six of these 23 SNVs had MAF <5%, and for these we performed additional follow-up in the GERA study. However, none reached genome-wide significance in a combined analysis of stage 1, EPIC-InterAct, and GERA (a total of 39,219 case and 198,990 control subjects) (Supplementary Table 6). Therefore, despite substantially enlarged sample sizes that would have allowed us to detect low-frequency risk alleles with modest effect sizes, the overwhelming majority of variants for which T2D association can be detected with these sample sizes are themselves common.

To identify loci containing multiple distinct signals, we performed approximate conditional analysis within the established and novel GWAS loci and detected two such novel common variant signals (Supplementary Table 7) (19,20). At the ANKRD55 locus, we identified a previously unreported distinct (Pconditional < 10−5) association signal led by rs173964 (Pconditional = 3.54 × 10−7, MAF = 26%) (Supplementary Table 7 and Supplementary Fig. 4). We also observed multiple signals of association at loci with previous reports of such signals (4,14), including CDKN2A/B (three signals in total), DGKB and KCNQ1 (six signals), and HNF4A and CCND2 (three signals) (Supplementary Table 7 and Supplementary Fig. 4). At CCND2, in addition to the main signal with lead SNV rs4238013, we detected 1) a novel distinct signal led by a common variant, rs11063018 (Pconditional = 2.70 × 10−7, MAF = 19%) and 2) a third distinct signal led by a low-frequency protective allele (rs188827514, MAF = 0.6%; ORconditional = 0.60, Pconditional = 1.24 × 10−6) (Supplementary Fig. 5A and Supplementary Table 7), which represents the same distinct signal as that at rs76895963 (Pconditional = 1.0) reported in the Icelandic population (Supplementary Fig. 5) (7). At HNF4A, we confirmed recent analyses (obtained in partially overlapping data) (14) that a low-frequency missense variant (rs1800961, p.Thr139Ile, MAF = 3.7%) is associated with T2D and is distinct from the known common variant GWAS signal (which we mapped here to rs12625671).

We evaluated the trans-ethnic heterogeneity of allelic effects (i.e., discordance in the direction and/or magnitude of estimated ORs) at novel loci on the basis of Cochran's Q statistics from the largest T2D trans-ancestry GWAS meta-analysis to date (2). Using reported summary statistics from that study, we observed no significant evidence of heterogeneity of effect size (Bonferroni correction PCochran's < 0.05/13 = 0.0038) between major ancestral groups at any of the 13 loci (Supplementary Table 8). These results are consistent with these loci being driven by common causal variants that are widely distributed across populations.

1000G Variant Density for Identification of Potentially Causal Genetic Variants

We used credible set fine-mapping (33) to investigate whether 1000G imputation allowed us to better resolve the specific variants driving 95 distinct T2D association signals at 82 loci (Supplementary Material). The 99% credible sets included between 1 and 7,636 SNVs; 25 included fewer than 20 SNVs, 16 fewer than 10 (Supplementary Tables 9 and 10). We compared 1000G-based credible sets with those constructed from HapMap SNVs alone (Fig. 2 and Supplementary Table 9). At all but three of the association signals (two at KCNQ1 and rs1800961 at HNF4A), 1000G imputation resulted in larger credible sets (median increase of 34 variants) spanning wider genomic intervals (median interval size increase of 5 kb) (Fig. 2 and Supplementary Table 9). The 1000G-defined credible sets included >85% of the SNVs in the corresponding HapMap sets (Supplementary Table 9). Despite the overall larger credible sets, we asked whether 1000G imputation enabled an increase in the posterior probability afforded to the lead SNVs, but we found no evidence to this effect (Fig. 2).

Figure 2

A: The number (N) of SNVs included in 99% credible sets when performed on all SNVs compared with when analyses were restricted to those SNVs present in HapMap. B: The cumulative πc of the top three SNVs among all 1000G SNVs and after restriction to HapMap SNVs is shown. While the low-frequency SNV at TP53INP1 (rs11786613) did not reach the threshold for a distinct signal in approximate conditional analyses, we fine-mapped both this variant and the previous common signal separately after reciprocal conditioning, which suggested they were independent. C: The MAF of the lead SNV identified in current analyses compared with that identified among SNVs present in HapMap. D: The association of the low-frequency variant rs11786613 (blue) and that of the previous lead variant at this locus, rs7845219 (purple). The low-frequency variant overlaps regulatory annotations active in pancreatic islets, among other tissues, and the sequence surrounding the A allele of this variant has an in silico recognition motif for a FOXA1:AR (androgen receptor) protein complex.

Within the 50 loci previously associated with T2D in Europeans (4), which had at least modest evidence of association in the current analyses (P < 5 × 10−4), we asked whether the lead SNV in 1000G-imputed analysis was of similar frequency to that observed in HapMap analyses. Only at TP53INP1 was the most strongly associated 1000G-imputed SNV (rs11786613, OR = 1.21, P = 1.6 × 10−6, MAF = 3.2%) of substantially lower frequency than the lead HapMap-imputed SNV (3) (rs7845219, MAF = 47.7%) (Fig. 2). rs11786613 was neither present in HapMap nor on the Metabochip (Supplementary Fig. 6). Reciprocal conditioning of this low-frequency SNV and the previously identified common lead SNV (rs7845219, OR = 1.05, P = 5.0 × 10−5, MAF = 47.5%) indicated that the two signals were likely to be distinct but the signal at rs11786613 did not meet our threshold (Pconditional < 10−5) for locus-wide significance (Supplementary Fig. 4).

Pathophysiological Insights From Novel T2D Associations

Among the 13 novel T2D-associated loci, many (such as those near HLA-DQA1, NRXN3, GIP, ABO, and CMIP) included variants previously implicated in predisposition to other diseases and traits (r2 > 0.6 with the lead SNV) (Supplementary Table 3 and Supplementary Material). For example, the novel association at SNV rs1182436 lies ∼120 kb upstream of MNX1, a gene implicated in pancreatic hypoplasia and neonatal diabetes (34–36).

The lead SNV rs78761021 at the GLP2R locus, encoding the receptor for glucagon-like peptide 2, is in strong LD (r2 = 0.87) with a common missense variant in GLP2R (rs17681684, D470N, P = 3 × 10−7). These signals were strongly dependent and mutually extinguished in reciprocal conditional analyses, consistent with the coding variant being causal and implicating GLP2R as the putative causal gene (Supplementary Fig. 7). While previously suggested to regulate energy balance and glucose tolerance (37), GLP2R has primarily been implicated in gastrointestinal function (38,39). In contrast, GLP1R, encoding the glucagon-like peptide 1 receptor (the target for a major class of T2D therapies [40]), is more directly implicated in pancreatic islet function, and variation at this gene has been associated with glucose levels and T2D risk (41).

We also observed associations with T2D centered on rs9271774 near HLA-DQA1 (Table 1), a region showing a particularly strong association with T1D (42). There is considerable heterogeneity within, and overlap between, the clinical presentations of T1D and T2D, but these can be partially resolved through measurement of islet cell autoantibodies (43). Such measures were not uniformly available across studies contributing to our meta-analysis (Supplementary Table 1). We therefore considered whether the adjacency between T1D and T2D risk loci was likely to reflect misclassification of individuals with autoimmune diabetes as case subjects in the current study.

Three lines of evidence make this unlikely. First, the lead T1D-associated SNV in the HLA region (rs6916742) was only weakly associated with T2D in the current study (P = 0.01), and conditioning on this variant had only modest impact on the T2D association signal at rs9271774 (Punconditional = 3.3 × 10−7; Pconditional = 9.1 × 10−6). Second, of 52 published genome-wide significant T1D association GWAS signals, 50 were included in the current analysis: only six of these reached even nominal association with T2D (P < 0.05; Supplementary Fig. 8), and at one of these six (BCAR1), the T1D risk allele was protective for T2D. Third, in genetic risk score analyses, the combined effect of these 50 T1D signals on T2D risk was of only nominal significance (OR = 1.02 [95% CI 1.00–1.03], P = 0.026), and significance was eliminated when the six overlapping loci were excluded (OR = 1.00 [95% CI 0.98–1.02], P = 0.73). In combination, these findings argue against substantial misclassification and indicate that the signal at HLA-DQA1 is likely to be a genuine T2D signal.

Potential Genes and Pathways Underlying the T2D Loci: eQTL and Pathway Analysis

cis-eQTLs analyses highlighted four genes as possible effector transcripts: ABO (pancreatic islets), PLEKHA1 (whole blood), and HSD17B12 (adipose, liver, muscle, whole blood) at the respective loci and HLA-DRB5 expression (adipose, pancreatic islets, whole blood) at the HLA-DQA1 locus (Supplementary Table 11).

We next asked whether large-scale gene expression data, mouse phenotypes, and protein–protein interaction networks could implicate specific gene candidates and gene sets in the etiology of T2D. Using DEPICT (32), 29 genes were prioritized as driving observed associations (FDR <0.05), including ACSL1 and CMIP among the genes mapping to the novel loci (Supplementary Table 12). These analyses also identified 20 enriched reconstituted gene sets (FDR <5%) falling into four groups (Supplementary Fig. 9) (complete results, including gene prioritization, can be downloaded from http://diagram-consortium.org/2017_Scott_DIAGRAM_1000G/). These included pathways related to mammalian target of rapamycin (mTOR) based on coregulation of the IDE, TLE1, SPRY2, CMIP, and MTMR3 genes (44).

Overlap of Associated Variants With Regulatory Annotations

We observed significant enrichment for T2D-associated credible set variants in pancreatic islet active enhancers and/or promoters (log odds [β] = 0.74, P = 4.2 × 10−8) and FOXA2 binding sites (β = 1.40, P = 4.1 × 10−7), as previously reported (Supplementary Table 13) (14). We also observed enrichment for T2D-associated variants in coding exons (β = 1.56, P = 7.9 × 10−5), in EZH2-binding sites across many tissues (β = 1.35, P = 5.3 × 10−6), and in binding sites for NKX2.2 (β = 1.73, P = 4.1 × 10−8) and PDX1 (β = 1.46, P = 7.4 × 10−6) in pancreatic islets (Supplementary Fig. 10).

Even though credible sets were generally larger, analyses performed on the 1000G imputed results produced stronger evidence of enrichment than equivalent analyses restricted to SNVs present in HapMap. This was most notably the case for variants within coding exons (β = 1.56, P = 7.9 × 10−5 in 1000G compared with β = 0.68, P = 0.62 in HapMap) and likely reflects more complete capture of the true causal variants in the more densely imputed credible sets. Single lead SNVs overlapping an enriched annotation accounted for the majority of the total posterior probability (πc > 0.5) at seven loci. For example, the lead SNV (rs8056814) at BCAR1 (πc = 0.57) overlaps an islet enhancer (Supplementary Fig. 11), while the newly identified low-frequency signal at TP53INP1 overlaps an islet promoter element (rs117866713, πc = 0.53) (Fig. 2) (31).

We applied hierarchical clustering to the results of diabetes-related quantitative trait associations for the set of T2D-associated loci from the current study, identifying three main clusters of association signals with differing impact on quantitative traits (Supplementary Table 9). The first, including GIPR, C2CDC4A, CDKAL1, GCK, TCF7L2, GLIS3, THADA, IGF2BP2, and DGKB, involved loci with a primary impact on insulin secretion and processing (26,29). The second cluster captured loci (including PPARG, KLF14, and IRS1) disrupting insulin action. The third cluster, showing marked associations with BMI and lipid levels, included NRXN3, CMIP, APOE, and MC4R but not FTO, which clustered alone.

In regulatory enhancement analyses, we observed strong tissue-specific enrichment patterns broadly consistent with the phenotypic characteristics of the physiologically stratified locus subsets. The cluster of loci disrupting insulin secretion showed the most marked enrichment for pancreatic islet regulatory elements (β = 0.91, P = 9.5 × 10−5). In contrast, the cluster of loci implicated in insulin action was enriched for annotations from adipocytes (β = 1.3, P = 2.7 × 10−11) and monocytes (β = 1.4, P = 1.4 × 10−12), and that characterized by associations with BMI and lipids showed preferential enrichment for hepatic annotations (β = 1.15, P = 5.8 × 10−4) (Fig. 3). For example, at the novel T2D-associated CMIP locus, previously associated with adiposity and lipid levels (28,45), the lead SNV (rs2925979, πc = 0.91) overlaps an active enhancer element in both liver and adipose tissue, among others (Supplementary Fig. 11).

Figure 3

T2D loci stratified by patterns of quantitative trait (e.g., glycemic, insulin, lipid, and anthropometric) effects show distinct cell-type annotation patterns. We hierarchically clustered loci based on endophenotype data and identified groups of T2D loci associated with measures of insulin secretion (A), insulin resistance (B), and BMI/lipids (C). We then tested the effect of variants in cell-type enhancer and promoter chromatin states on the posterior probabilities of credible sets for each group. We identified most significant effects among pancreatic (Panc.) islet chromatin for insulin secretion loci, CD14+ monocyte and adipose chromatin for insulin resistance loci, and liver chromatin for BMI/lipid loci.

Discussion

In this large-scale study of T2D genetics, in which individual variants were assayed in up to 238,209 subjects, we identified 13 novel T2D-associated loci at genome-wide significance and refined causal variant location for the 13 novel and 69 established T2D loci. We also found evidence for enrichment in regulatory elements at associated loci in tissues relevant for T2D and demonstrated tissue-specific enrichment in regulatory annotations when T2D loci were stratified according to inferred physiological mechanism.

We calculate that the present analysis, together with loci reported in other recent publications (9), brings the total number of independent T2D associations to 128 distinct signals at 113 loci (Supplementary Table 3). Lead SNVs at all 13 novel loci were common (MAF >15%) and of comparable effect size (1.07 ≤ OR ≤ 1.10) to previously identified common variant associations (2,4). Associations at the novel loci showed homogeneous effects across diverse ethnicities, supporting the evidence for coincident common risk alleles across ancestry groups (2). Moreover, we conclude that misclassification of diabetes subtype is not a major concern for these analyses and that the HLA-DQA1 signal represents genuine association with T2D, independent of nearby signals that influence T1D.

We observed a general increase in the size of credible sets with 1000G imputation compared with HapMap imputation. This is likely due to improved enumeration of potential causal common variants on known risk haplotypes rather than resolution toward low-frequency variants of larger effect driving common variant associations. These findings are consistent with the inference (arising also from the other analyses reported here) that the T2D risk signals identified by GWAS are overwhelmingly driven by common causal variants. In such a setting, imputation with denser reference panels, at least in ethnically restricted samples, provides more complete elaboration of the allelic content of common risk haplotypes. Finer resolution of those haplotypes that would provide greater confidence in the location of causal variants will likely require further expansion of trans-ethnic fine-mapping efforts (2). The distinct signals at the established CCND2 and TP53INP1 loci point to contributions of low-frequency variant associations of modest effect but indicate that even larger samples will be required to robustly detect association signals at low frequency. Such new large data sets might be used to expand the follow-up of suggestive signals from our analysis.

The discovery of novel genome-wide significant association signals in the current analysis is attributable primarily to increased sample size rather than improved genomic coverage. Although we queried a large proportion of the low-frequency variants present in the 1000G European reference haplotypes and had >80% power to detect genome-wide significant associations with OR >1.8 for the tested low-frequency risk variants, we found no such low-frequency variant associations in either established or novel loci. While low-frequency variant coverage in the current study was not complete, this observation adds to the growing evidence (2,4,9,46) that few low-frequency T2D risk variants with moderate to strong effect sizes exist in European ancestry samples and is consistent with a primary role for common variants of modest effect in T2D risk. The current study reinforces the conclusions from a recent study that imputed from whole-genome sequencing data—from 2,657 European T2D case and control subjects rather than 1000G—into a set of GWAS studies partially overlapping with the present meta-analysis. We demonstrated that the failure to detect low-frequency associations in that study is not overcome by a substantial increase in sample size (9). It is worth emphasizing that we did not, in this study, have sufficient imputation quality to test for T2D associations with rare variants and we cannot evaluate the collective contribution of variants with MAF <0.5% to T2D risk.

The development of T2D involves dysfunction of multiple mechanisms across several distinct tissues (9,29,31,47,48). When coupled with functional data, we saw larger effect estimates for enrichment of coding variants than observed with HapMap SNVs alone, consistent with more complete recovery of the causal variants through imputation using a denser reference panel. The functional annotation analyses also demonstrated that the stratification of T2D risk loci according to primary physiological mechanism resulted in evidence for consistent and appropriate tissue-specific effects on transcriptional regulation. These analyses exemplify the use of a combination of human physiology and genomic annotation to position T2D GWAS loci with respect to the cardinal mechanistic components of T2D development. Extension of this approach is likely to provide a valuable in silico strategy to aid prioritization of tissues for mechanistic characterization of genetic associations. Using the hypothesis-free pathway analysis of T2D associations with DEPICT (32), we highlighted a causal role of mTOR signaling pathway in the etiology of T2D not observed from individual loci associations. The mTOR pathway has previously been implicated in the link between obesity, insulin resistance, and T2D from cell and animal models (44,49).

The current results emphasize that progressively larger sample sizes, coupled with higher density sequence-based imputation (13), will continue to represent a powerful strategy for genetic discovery in T2D and in complex diseases and traits more generally. At known T2D-associated loci, identification of the most plausible T2D causal variants will likely require large-scale multiethnic analyses, where more diverse haplotypes, reflecting different patterns of LD, in combination with functional (31,50,51) data allow refinement of association signals to smaller numbers of variants (2).