Identification of type 2 diabetes loci in 433,540 East Asian individuals.

Meta-analyses of genome-wide association studies (GWAS) have identified more than 240 loci that are associated with type 2 diabetes (T2D)1,2; however, most of these loci have been identified in analyses of individuals with European ancestry. Here, to examine T2D risk in East Asian individuals, we carried out a meta-analysis of GWAS data from 77,418 individuals with T2D and 356,122 healthy control individuals. In the main analysis, we identified 301 distinct association signals at 183 loci, and across T2D association models with and without consideration of body mass index and sex, we identified 61 loci that are newly implicated in predisposition to T2D. Common variants associated with T2D in both East Asian and European populations exhibited strongly correlated effect sizes. Previously undescribed associations include signals in or near GDAP1, PTF1A, SIX3, ALDH2, a microRNA cluster, and genes that affect the differentiation of muscle and adipose cells3. At another locus, expression quantitative trait loci at two overlapping T2D signals affect two genes-NKX6-3 and ANK1-in different tissues4-6. Association studies in diverse populations identify additional loci and elucidate disease-associated genes, biology, and pathways.

Type 2 diabetes (T2D) is a common metabolic disease primarily caused by insufficient insulin production and/or secretion by the pancreatic β cells and insulin resistance in peripheral tissues[7]. Most genetic loci associated with T2D have been identified in populations of European (EUR) ancestry, including a recent meta-analysis of genome-wide association studies (GWAS) of nearly 900,000 individuals of European ancestry that identified >240 loci influencing T2D risk[1]. Differences in allele frequency between ancestries affect the power to detect associations within a population, particularly among variants rare or monomorphic in one population but more frequent in another[2,8]. Although smaller than studies in European populations, a recent T2D meta-analysis in almost 200,000 Japanese individuals identified 28 additional loci[2]. The relative contributions of different pathways to T2D pathophysiology may also differ between ancestry groups. For example, in East Asian (EAS) populations, T2D prevalence is greater than in European populations among people of similar body mass index (BMI) or waist circumference[9]. To identify new genetic associations and provide insight into T2D pathogenesis, we performed the largest meta-analysis of East Asian individuals to date.

RESULTS

We conducted a fixed-effect inverse-variance weighted GWAS meta-analysis combining 23 studies imputed to the 1000 Genomes Phase 3 reference panel from the Asian Genetic Epidemiology Network (AGEN) consortium (Supplementary Tables 1–3). We performed sex-combined T2D association without BMI adjustment in 77,418 T2D cases and 356,122 controls (effective sample size, Neff=211,793). For the subset of 54,481 T2D cases and 224,231 controls (Neff= 135,780) with BMI available, additional analyses were performed with and without BMI adjustment in sex-combined and sex-stratified models (Extended Data Figure 1). We defined “lead” variants as the strongest T2D-associated variants with P<5x10−8 and defined the region +/− 500 kb from the lead variant as a locus. A locus was considered novel if the lead variant was located at least 500 kb from previously reported T2D-associated variants in any ancestry.

Extended Data Figure 1:

Flow chart of study design, depicting the different data analyses performed.

Using summary association statistics for ~11.8 million variants, without adjustment for BMI (Extended Data Figure 1; Supplementary Tables 1–3), we identified lead variants associated with T2D at 183 loci, of which 51 were novel (Extended Data Table 1; Extended Data Figure 2; Supplementary Table 4). Lead variants at all novel loci were common (MAF≥5%; Extended Data Figure 3), except for low-frequency variants near GDAP1 (MAF=2.4%), which regulates mitochondrial proteins and metabolic flux in skeletal muscle[10], and PTF1A (MAF=4.7%), which encodes a transcription factor required for pancreatic acinar cell development[11]. Lead variants met a stricter P-value threshold for significance based on Bonferroni correction for 11.8 million tests (P<4.2x10−9) at 146 of the 183 loci, including 29 of the 51 novel loci.

Extended Data Table 1:

Novel lead variants associated with type 2 diabetes in East Asians.

Single variant association results from East Asian fixed-effect inverse-variance meta-analysis (BMI-unadjusted sex-combined model adjusted for age, sex, and study-specific covariates) using METAL. Loci were defined as novel if the lead variant is (1) at least 500 kb away and confirmed by GCTA to be conditionally independent from previously reported T2D-associated variants in any ancestry, and (2) assessed using locuszoom plots and biology lookups to be to be away from known loci with extended LD. Four additional variants met the definition for a novel locus but are located within the previously reported major histocompatibility complex (MHC) region; see Supplementary Table S5 for the full list of distinct association signals at the MHC region. rs4804181 was >500kb from primary signal rs3111316 in European meta-analysis, but <500kb from their secondary signal rs755734872. Genome-wide significant association is defined as P < 5.0E-08. Physical position based on hg19. Effect alleles are associated with increased risk for T2D. Odds ratios reflect per allele effects of variants on T2D risk.

Chr, chromosome; Pos, position; RAF, risk allele frequency; Neff, effective sample size; OR, odds ratio; CI, confidence interval; P, P-value

Extended Data Figure 2:

Manhattan plot for East Asian T2D meta-analysis association results in model unadjusted for BMI.

−log10(P) values from two-sided fixed-effects inverse-variance genome-wide meta-analysis association results for each variant (y-axis; maximal Neff=211,793) was plotted against the genomic position (hg19; x-axis). Known T2D loci achieving genome-wide significance (P<5.0x10−8) meta-analysis are shown in blue. Loci achieving genome-wide significance that are previously unreported for T2D association are shown in red.

Extended Data Figure 3:

The relationship between effect size and minor allele frequency.

Odds ratios (y-axis) and minor allele frequencies (x-axis) for 189 primary association signals from the T2D BMI-unadjusted models. Odds ratios are from two-sided fixed-effects inverse-variance meta-analysis on a maximal effective sample size of 211,793.

Using GCTA[12], we identified 301 distinct association signals that met a locus-wide significance threshold of P<1x10−5 (Supplementary Table 5), 228 of which were genome-wide significant (P<5x10−8). Overall, we observed 2-4 signals at 46 loci and ≥5 signals at 12 loci. Among the ten loci with the most significant meta-analysis P-values of association, seven contained ≥5 distinct signals (17 signals at INS/IGF2/KCNQ1; 7 signals at CDKN2A/B and GRM8/PAX4/LEP; 5 signals at CDKAL1, HHEX/IDE, CDC123/CAMK1D, and TCF7L2; Extended Data Figure 4; Supplementary Table 5). The seven signals at the GRM8/PAX4/LEP locus span 1.4 Mb, and no evidence of T2D association at this locus has yet been reported in non-East Asian ancestry groups[1,13] (Extended Data Figure 4C). Joint analyses confirmed independent associations (LD r2=0.0025) at two previously reported PAX4 missense variants[14], rs2233580 [Arg192His: risk allele frequency (RAF)=8.6%, OR=1.31, 95% CI 1.28 – 1.34, PGCTA=3.4x10−93] and rs3824004 (Arg192Ser: RAF=3.4%, OR=1.24, 95% CI 1.19-1.28, PGCTA=1.1x10−30). The association signals at this locus also include variants near LEP, which encodes leptin, a hormone that regulates appetite[15]; increased leptin levels are associated with obesity and T2D[16].

Extended Data Figure 4:

Regional association plots at three T2D-associated loci with the strongest association P-value and more than five distinct association signals in East Asians.

(A) INS/IGF2/KCNQ1, (B) CDKN2A/B, (C) PAX4/LEP. −log10(P) values were from the two-sided fixed-effect inverse-variance meta-analysis. Distinct signals (P<1.0x10−6 from GCTA conditional analyses) were plotted; Neff for each distinct signal are reported in Supplementary Table 4. Variants are colored based on East Asian 1000G Phase 3 LD with the lead variants for each association signal, shown as diamonds.

At the previously reported ANK1/NKX6-3 locus[1,17], we observed three distinct T2D association signals, two of which overlap and consist of variants spanning only ~25 kb (Figure 1). Given conflicting interpretation of candidate genes[1,5,18], we compared the T2D-association signals identified in East Asian individuals to eQTLs reported at this locus in islets[1,18-20], subcutaneous adipose[6], and skeletal muscle[5]. At the strongest signal, the lead T2D-associated variant rs33981001 is in high LD with the lead cis-eQTL variant for NKX6-3 in pancreatic islets (rs12549902; EAS LD r2=0.79, EUR r2=0.83)[18], and the T2D risk allele is associated with decreased expression of NKX6-3 (β=−0.36, P=6.1x10−7; Figure 1)[4]. NKX6-3, or NK6 homeobox 3, encodes a pancreatic islet transcription factor required for the development of alpha and β cells in the pancreas[21] and has been shown to influence insulin secretion[16]. At the second T2D-association signal, rs62508166 is in high LD with the lead cis-eQTL variant for ANK1 in subcutaneous adipose tissue[19] and skeletal muscle[15] (rs516946; EAS LD r2=0.96, EUR r2=0.80), and the T2D risk allele is associated with increased expression of ANK1 (subcutaneous adipose: β=0.20, P=1.8x10−7; skeletal muscle: β=1.01, P=2.8x10−22). ANK1 belongs to the ankyrin family of integral membrane proteins and has been shown to affect glucose uptake in skeletal muscle, and changes in its expression level may lead to insulin resistance[22]. Together, these GWAS and eQTL signals suggest that variants within this ~25 kb region act to increase or decrease expression levels of two different genes in different tissues to increase T2D risk.

Figure 1:

Two distinct T2D-association signals at the ANK1-NKX6-3 locus associated with expression levels of two transcripts in two tissues.

(A) Regional association plot for East Asian sex-combined BMI-unadjusted two-sided fixed-effect inverse-variance meta-analysis at ANK1-NKX6-3 locus. Approximate conditional analysis using GTCA identified three distinct T2D-association signals (P<1x10−5) at this locus (signal 1, rs33981001, Neff=211,793; signal 2, rs62508166, Neff=211,793; signal 3, rs144239281, Neff=208,431, in order of strength of association). Using 1000G Phase3 East Asian LD, variants are colored in red and blue with the first and second distinct signals respectively (lead variants represented as diamonds). (B) Variant rs12549902, in high LD (EAS LD r2=0.80, EUR r2=0.83) with T2D signal 1, shows the strongest association with expression levels of NXK6-3 in pancreatic islets in 118 individuals[16]. (C) Variant rs516946, in high LD (EAS LD r2=0.96, EUR r2=0.80) with T2D signal 2, shows the strongest association with expression levels of ANK1 in subcutaneous adipose tissue in 770 individuals[19]. As rs62508166 is not available in the subcutaneous adipose tissue data set, a variant in perfect LD (rs28591316) was used and is represented by the blue diamond.

In T2D association analyses adjusted for BMI, we identified an additional six loci, four of which were not reported previously for T2D, including loci near MYOM3/SRSF10, TSN, GRB10, and NID2 (Supplementary Table 4). At the NID2 locus, the T2D risk allele is very rare or monomorphic in non-East Asian individuals and has previously demonstrated significant associations with lower BMI and higher triglycerides in East Asian individuals, consistent with a lipodystrophy phenotype[23,24]. The lead GRB10 variant is in low LD (EUR r=0.08, EAS r=0.57) with a variant associated with glucose-stimulated insulin secretion in European individuals[25].

Across the models with and without adjustment for BMI, correlation for the effect sizes genome-wide was higher in East Asian individuals (r=0.98) than in European individuals (r=0.89). For the 189 T2D-associated loci in East Asian individuals, the correlation increased to 0.99 (Extended Data Figure 5). Loci with larger effects in BMI-adjusted models include FGFR2 and NID2, identified only in East Asian populations and associated with lipodystrophy traits or body fat distribution. These results may reflect the role of body fat distribution in insulin resistance and T2D among East Asian individuals.

Extended Data Figure 5:

Effect size comparison of lead variants in sex-combined models unadjusted and adjusted for BMI.

At 189 lead variants identified in the East Asian BMI-unadjusted sex-combined T2D meta-analysis, per-allele effect sizes (β) from the BMI-adjusted sex-combined model were plotted against the BMI-unadjusted sex-combined model. Both sex-combined models were from two-sided fixed-effect inverse-variance meta-analyses and included the same set of studies for comparable sample size. Each point denotes the per-allele effect size; standard errors of the effect size estimates extend out as grey lines. Effect sizes between the two models are highly correlated with a Pearson correlation coefficient r=0.99 (Supplementary Table 4).

In sex-stratified analyses of males (28,027 cases and 89,312 controls) and females (27,370 cases and 135,055 controls), we identified six additional novel sex-specific loci: (i) three male-specific loci near FOXK1, PDE3A, and IFT81, and one female-specific locus near LMTK2 in models without adjustment for BMI, and (ii) one male-specific locus near LINC00851 and one female-specific locus near CPS1 in models with adjustment for BMI (Supplementary Table 6). The lead CPS1 variant rs1047891 (Thr1412Asn) has been reported to have a stronger effect in females than in males for cardiovascular disease and several blood metabolites[26]. Taken together, we identified a total of 61 novel loci across BMI-unadjusted, BMI-adjusted, and sex-stratified models, of which 33 met a stricter P-value threshold (P<4.2x10−9).

Among all T2D-associated loci, a region spanning ~2 Mb near ALDH2 exhibited the strongest differences between sexes (rs12231737, Phet=2.6x10−19), with compelling evidence of association in males (Pmales=5.8x10−27) and no evidence for association in females (Pfemales=0.19) (Extended Data Figure 6; Supplementary Table 6). This sex difference is also observed after adjusting for BMI (Pmales_adjBMI=5.2x10−21, Pfemales_adjBMI=0.053). Further, joint conditional analyses revealed two conditionally distinct signals (rs12231737, PGCTA=1.7x10−21; rs557597782, PGCTA=4.9x10−7) in males only. ALDH2 encodes aldehyde dehydrogenase 2 family member, a key enzyme in alcohol metabolism that converts acetaldehyde into acetic acid. This stretch of T2D associations in males reflects a long LD block that arose due to a recent selective sweep in East Asian individuals and results in flushing, nausea, and headache following alcohol consumption[27]. The most significantly associated missense variant in moderate LD with rs12231737 (r2=0.68) was rs671 (ALDH2 Glu504Lys: RAF=77.7%, OR=1.17, 95% CI 1.14 – 1.20, Pmales=1.5x10−24), which leads to reduced ALDH2 activity and reduced alcohol metabolism, and has been associated with cardiometabolic traits in East Asian populations. The T2D risk allele is associated with better tolerance for alcohol; increased BMI, blood pressure, and high-density lipoprotein cholesterol; and decreased low-density lipoprotein cholesterol and cardiovascular risk[28-30]. The strong sexual dimorphism observed at this locus may be due to differences in alcohol consumption patterns between males and females[28,30], effects of alcohol on BMI, and/or differences in the effect of alcohol on insulin sensitivity[31].

Extended Data Figure 6:

Regional plots of male-specific T2D-associated locus, ALDH2.

For each plot, −log10(P) values from association results from two-sided fixed-effect inverse-variance meta-analyses for each variant (y-axis) was plotted against the genomic position (hg19; x-axis). The lead variant rs12231737 plotted is the lead variant from the BMI-unadjusted male-specific meta-analysis (Neff=65,202) and also the sex-combined meta-analysis (Neff=138,947) from the same subset of individuals included in the sex-stratified analyses (female-specific Neff=70,051). This lead variant rs12231737 is in high LD with rs77768175, identified from the larger BMI-unadjusted sex-combined meta-analysis (East Asian r2=0.80). (A) Males only, (B) sex-combined, and (C) females only. Variants are shaded based on East Asian 1000G Phase 3 LD with the lead variant, shown as a purple diamond.

With an effective sample size comparable to the largest study of T2D in European individuals (East Asian Neff=211,793; European Neff = 231,436)[1] and imputation to a dense 1000 Genomes reference panel, our results provide the most comprehensive and precise catalogue of East Asian T2D effects to date for comparisons across ancestries (Figure 2; Supplementary Table 7). For 183 EAS T2D loci and 231 EUR T2D loci (unadjusted for BMI)[1], we compared the per-allele effect sizes for the 332 variants available in both datasets (i.e. polymorphic and passed quality control), including lead variants from both ancestries at shared signals. Overall, the per-allele effect sizes between the two ancestries were moderately correlated (r=0.55; Figure 2A). When the comparison was restricted to the 278 variants that are common (MAF≥5%) in both ancestries, the effect size correlation increased to r=0.59 (Figure 2B; Extended Data Figure 7). This effect size correlation further increased to r=0.87 for 106 variants significantly associated with T2D (P<5x10−8) in both ancestries. Based on Cochran’s heterogeneity test, 28 of 332 variants (8.4%) exhibited significant heterogeneity in effect sizes between East Asian and European populations, including 22 that were significant in only one population (Supplementary Table 7) and six with larger effect sizes in one population (e.g. CDKAL1, KCNQ1, and HNF1B). While the overall effect sizes for all 332 variants appear, on average, to be stronger in East Asian individuals than European individuals, this trend is reduced when each locus is represented only by the lead variant from one population (Extended Data Figure 8). Specifically, 39 variants identified in the European meta-analysis with imputation using the Haplotype Reference Consortium panel are missing from the comparison because they were rare/monomorphic or poorly imputed in the East Asian meta-analysis, with imputation based on the smaller and more heterogenous 1000 Genomes reference panel.

Figure 2:

Effect size comparison of lead variants identified in this East Asian T2D GWAS BMI-unadjusted meta-analysis and previous European T2D GWAS meta-analysis[2].

For 332 unique lead variants identified from the two BMI unadjusted meta-analyses, per-allele effect sizes (β) from the European meta-analysis (y-axis) were plotted against per-allele effect sizes from this East Asian meta-analysis (x-axis). Effect sizes from both meta-analyses were from two-sided fixed-effect inverse-variance meta-analysis (maximal Neff=211,793 for East Asian and 231,436 for European meta-analyses). Each point denotes the per-allele effect size; standard errors of the effect size estimates extend out as grey lines. (A) All 332 lead variants; (B) 278 lead variants with minor allele frequency ≥5% in both ancestries. Variants are colored purple if they were significant (P<5x10−8) in the East Asian analysis only, green if they were significant in European analysis only, and blue if they were significant in both the East Asian and European analyses (see Methods and Supplementary Table 7). The dashed diagonal line represents the trend line across all plotted variants. Compared to Supplementary Table 7, 70 variants are not plotted; 31 variants were present only in the analysis of East Asian individuals (median effect size 0.065; interquartile range 0.049-0.110) and 39 variants were present only in the analysis of European individuals (median effect size 0.083; interquartile range 0.063-0.170).

Extended Data Figure 7:

Effect size comparison of common lead variants (MAF≥5%) identified in this East Asian meta-analysis and a previously published European T2D GWAS meta-analysis[2].

For 278 unique lead variants with MAF≥5% in both the East Asian and European BMI-unadjusted meta-analyses, per-allele effect sizes (β) from Mahajan et al.[2] (y-axis) were plotted against per-allele effect sizes from this East Asian meta-analysis (x-axis). Effect sizes from both meta-analyses were from two-sided fixed-effect inverse-variance meta-analyses (maximal Neff=211,793 for East Asian and 231,436 for European meta-analyses). Each point denotes the per-allele effect size; standard errors of the effect size estimates extend out as grey lines. Variants are colored purple if they were significant in the East Asian meta-analysis only, green if they were significant in European meta-analysis only, and blue if they were significant in both the East Asian and European meta-analyses. (see Methods and Supplementary Table 7).

Extended Data Figure 8:

Effect size comparison of lead variants identified in East Asian BMI-unadjusted meta-analysis and previously published European T2D GWAS meta-analysis[2].

For 332 lead variants identified from the two BMI-unadjusted meta-analyses, per-allele effect sizes (β) from a European meta-analysis (y-axis) were plotted against per-allele effect sizes from this East Asian meta-analysis (x-axis). Effect sizes from both meta-analyses were from two-sided fixed-effect inverse-variance meta-analysis (maximal Neff=211,793 for East Asian and 231,436 for European meta-analyses). Each point denotes the per-allele effect size; standard errors of the effect size estimates extend out as grey lines. (A) 152 lead variants significant in the East Asian meta-analysis (purple) or both the East Asian and European meta-analysis (blue) and (B) 192 lead variants significant in the European meta-analysis (green) or both the East Asian and European meta-analysis (blue). These plots include only one variant per locus, in contrast to Figure 2 and Extended Data Figure 7.

Variants exhibiting the largest differences in effect sizes across ancestries are generally rare (MAF ≤0.1%) in European populations but common (e.g. PAX4, RANBP3L) or low-frequency (e.g. ZNF257, DGKD) in East Asian populations. For example, rs142395395 near ZNF257 (RAF=96.9%, OR=1.24, 95% CI 1.19-1.29, P=7.0x10−23) has been reported only twice in 15,414 individuals of non-Finnish European ancestry from the gnomAD database[32]. This variant tags a previously described inversion of 415 kb observed only in East Asian individuals that disrupts the coding sequence and expression of ZNF257, as well as lymphoblastoid expression of 81 downstream genes and transcripts[33]. These data suggest that ZNF257 and/or downstream target genes influence T2D susceptibility.

We identified many loci for which the lead variants exhibited similar allele frequencies and effect sizes in both the East Asian and European meta-analyses, but only reached genome-wide significance in the East Asian meta-analysis. Given shared susceptibility across ancestry groups, these loci may be detected in non-East Asian populations when sample sizes increase. Among these variants is rs117624659, located near NKX6-1 (PEAS = 2.0x10−16, PEUR =2.2x10−4). This lead variant overlaps a highly conserved region that shows open chromatin specific to pancreatic islets. We conducted transcriptional reporter assays in MIN6 mouse insulinoma cells and observed that rs117624659 exhibited significant allelic differences in enhancer activity (Figure 3). In the pancreas, NK6 homeobox 1 (NKX6.1) is required for the development of insulin-producing β cells and is a potent bifunctional transcriptional regulator[34]. Further, inactivation of Nkx6.1 in mice demonstrated rapid-onset diabetes due to defects in insulin biosynthesis and secretion[35]. Unexpectedly, the T2D risk allele showed increased transcriptional activity, suggesting that the variant does not act in isolation or that NXK6-1 is not the target gene.

Figure 3:

rs117624659 at NKX6-1 locus exhibits allelic differences in transcriptional activity.

(A) rs117624659 (Neff=211,214; purple diamond) shows the strongest association with T2D in the region. P values were from two-sided fixed-effect inverse-variance meta-analysis. Variants are colored based on 1000G Phase 3 East Asian LD with rs117624659. (B) rs117624659 and an additional candidate variant rs142390274 in high pairwise LD (r2>0.80) span a 22 kb region approximately 75 kb upstream of NKX6-1. rs117624659 overlaps a region of open chromatin in pancreatic islets and lies within a region conserved across vertebrates. (C) rs117624659-T, associated with increased risk of T2D, showed greater transcriptional activity in an element cloned in both forward and reverse orientations with respect to NKX6-1 in MIN6 cells compared to rs117624659-C and an “empty vector” containing a minimal promoter. Black lines represent mean (center horizontal line) and standard error (extended lines) relative luciferase activity from two-sided, unpaired t-tests using data from n=5 biologically independent samples/independent experiments.

At one of the novel T2D-associated loci near SIX3, the risk allele of East Asian lead variant rs12712928-C (RAF=40.2%, OR=1.06, 95% CI 1.04 – 1.07, P=1.8x10−14) is common across non-East Asian ancestries, ranging from 16.0% in Europeans to 26.4% in South Asians; however, there was no evidence of association in the other ancestry groups (meta-analysis: OR=0.98, 95% CI 0.96 – 0.99, P=2.9x10−3; Extended Data Figure 9A, Supplementary Table 8). Within the East Asian meta-analysis, the direction of effect is consistent across East Asian countries (Extended Data Figure 9B) and within the contributing cohorts (Extended Data Figure 9C). The T2D risk allele rs12712928-C is associated with higher fasting glucose levels in East Asian populations, has the strongest association with lower expression levels of both SIX3 and SIX2 in pancreatic islets[19], and demonstrated allele-specific binding to the transcription factor GABPA and significantly lower levels of transcriptional activity[36]. While rs12712928-C is present on only one common haplotype in most populations, it is present on an additional common haplotype (frequency =0.075) in East Asians, suggesting that the effect size attributed to rs12712928 may be influenced by other nearby unknown variants.

Extended Data Figure 9:

Forest plots of BMI-unadjusted meta-analysis association results at SIX3-SIX2 locus.

Odds ratios (black boxes) and 95% confidence intervals (horizontal lines extending out) for T2D associations at the lead East Asian variant (rs12712928) are presented (A) across ancestries of African-American (AFR), East Asian (EAS), European (EUR)[2], Hispanic (HIS), and South Asian (SAS) individuals, (B) within four major East Asian populations (Chinese, Japanese, Korean, and Malay/Filipino combined due to small sample sizes), (C) from each contributing cohort. Effect sizes from East Asian study, ancestry, population, and combined meta-analysis were from two-sided fixed-effect inverse-variance meta-analysis. The size of the box is proportional to the sample size of each contributing study/ancestry/population, which are available in Supplementary Table 8. This East Asian study had >90% power to detect the observed association with a MAF=0.40, OR=1.06, and 77,418 T2D cases. Given the number of T2D cases and frequency of rs12712928-C within the other datasets, at 80% power, we can reasonably exclude association OR >1.07 in EUR and >1.15 in AFR, HIS, and SAS between rs12782928 and T2D. Full study names can be found in Supplementary Table 1 and corresponding sample sizes can be found in Supplementary Table 2.

To identify potential candidate genes underlying the T2D-association signals identified in East Asian individuals, we further characterized 92 known and novel loci for which the lead variant at the primary East Asian association signal is located >500 kb from the lead variant of any European T2D association signal[2] (Supplementary Table 9). We characterized loci using prior trait associations, cis-regulatory effects on expression (colocalized eQTL), predicted effects on protein sequence, and a literature search (Supplementary Tables 10–13). Based on association results from cardiometabolic trait consortia[37], Biobank Japan[38], and the UK Biobank[39], the lead T2D-associated variant at 18 of the 88 loci was associated (P<5x10−8) with at least one additional cardiometabolic trait, most frequently BMI or a fat mass trait (15 loci; Supplementary Tables 10 and 12). At 12 of the examined loci, T2D signals were colocalized with cis-eQTLs for transcripts in subcutaneous adipose tissue (n=5), skeletal muscle (n=3), pancreas (n=2), islets (n=3), or blood (n=5; Supplementary Tables 11–12), generating hypotheses of target genes and directions of effect; further examination of these candidate genes is warranted. At 19 loci, the lead T2D-associated variant or a proxy (East Asian r2>0.80) alter the protein sequence (Supplementary Tables 12). These variants affect mesenchymal stem cell differentiation and adipogenesis (GIT2, STEAP2 and JMJD1C), muscle stem cell biology (CALCR), glucose metabolism (PGM1 and SCTR), and insulin secretion (FGFR4; Supplementary Table 13). At SCTR, which encodes the G-protein coupled secretin receptor, the lead variant encodes Ala122Pro, located in the hormone receptor domain. While mechanistic inference is required, these potential molecular mechanisms suggest new T2D susceptibility genes primarily detected by analyses in East Asian individuals.

T2D loci were also identified at clusters of noncoding RNAs with roles in islet β cell function. One locus includes a set of microRNAs specifically expressed in islet β cells, the maternally expressed noncoding RNA MEG3, and the paternally expressed gene DLK1. Targets of these microRNAs increase β cell apoptosis[40], and reduced Meg3 impairs insulin secretion[41]. DLK1 inhibits adipocyte differentiation, protecting from obesity[3], and promotes pancreatic ductal cell differentiation into β cells, increasing insulin secretion[42,43]. Other variants near MEG3 have been associated with type 1 diabetes (EAS and EUR LD r2=0 with EAS lead variant)[44]. The other noncoding RNA locus is the MIR17HG cluster of miRNAs that regulate glucose-stimulated insulin secretion and pancreatic β cell proliferation stress[45]; one of these microRNAs, miR-19a, affects hepatic gluconeogenesis[46]. Yet another T2D locus is located near TRAF3, which is a direct target of the MIR17HG microRNA cluster and promotes hyperglycemia by increasing hepatic glucose production[47,48]. The T2D association results suggest that these noncoding RNAs influence disease susceptibility.

DISCUSSION

These T2D GWAS meta-analyses in the largest number of East Asian individuals analyzed to date identified 61 novel loci, providing additional insight into the biological basis of T2D. The results emphasize substantial shared T2D susceptibility with European individuals, as shown by the strong correlation of effect sizes among T2D-associated genetic variants with common allele frequencies in both East Asian and European ancestry populations. Compared to a recent T2D study in individuals of European ancestry[1], we observed less attenuation of effects on T2D in analyses adjusted for BMI. Loci with a greater effect on T2D after adjusting for BMI include loci with lipodystrophy-like traits identified only in East Asian individuals to date, adding support to the observation[49,50] that factors beyond overall BMI, such as visceral adiposity or lipodystrophy, may also play a role in T2D in East Asians. The results also detect novel associations in East Asian individuals identified because they have higher allele frequencies in East Asian populations, exhibit larger effect sizes, and/or are influenced by other environmental or behavioral factors such as alcohol consumption.

The identified loci point to multiple plausible molecular mechanisms and many new candidate genes linking T2D susceptibility to diverse biological processes. Following the annotation of loci identified in the East Asian meta-analysis, we speculate a substantial role for insulin resistance in T2D pathogenesis among East Asian individuals through skeletal muscle, adipose, and liver development and function. We also provide evidence that multiple distinct association signals in the same region do not necessarily act through the same gene. Conditionally distinct association signals in close proximity can affect different genes that may act in different tissues by different mechanisms, emphasizing the value of identifying functional variants that enable variant-to-gene links to be examined directly. Our results provide a foundation for future biological research in T2D pathogenesis and offer potential targets for mechanisms for interventions in disease risk.

METHODS

Ethics statement

All human research was approved by the relevant institutional review boards for each study at their respective sites (Supplementary Table 1) and conducted according to the Declaration of Helsinki. All participants provided written informed consent.

Study cohorts and quality control

The East Asian type 2 diabetes (T2D) meta-analyses were performed with studies participating in the Asian Genetic Epidemiology Network (AGEN), a consortium of genetic epidemiology studies of T2D and related traits conducted in individuals of East Asian ancestry, and the Diabetes Meta-analysis of Trans-ethnic Association Studies (DIAMANTE), a consortium examining the genetic contribution to T2D across diverse ancestry populations including African-American, East Asian, European, Hispanic, and South Asian. The East Asian meta-analysis included 77,418 T2D cases and 356,122 controls from 23 GWAS, including three biobanks, CKB, KBA[51,52], and BBJ[2] [effective sample size (Neff) = 211,793; Extended Data Figure 1]. A subset of studies with BMI measurement available was analyzed with and without BMI adjustment in sex-combined and sex-specific models (54,481 cases, 224,231 controls; Neff = 135,780). For each study, T2D case control ascertainment is described in Supplementary Table 1 and summary statistics are provided in Supplementary Table 2. As T2D case definitions across cohorts differ, it is possible that cases of type 1 diabetes and maturity onset diabetes of the young (MODY) are included in these meta-analyses. Included studies were genotyped on either commercially available or customized Affymetrix or Illumina genome-wide genotyping arrays. Array quality control criteria implemented within each study, including variant call rate and Hardy-Weinberg equilibrium, are summarized in Supplementary Table 3. To harmonize study-level genotype scaffold for imputation to 1000 Genomes (1000G) reference panels, each study adopted a uniform protocol for pre-imputation quality checks. Each study applied the protocol to exclude variants with: i) mismatched chromosomal positions or alleles not present in the reference panel; ii) ambiguous alleles (AT/CG) with minor allele frequency (MAF) >40% in the reference panel; or iii) absolute allele frequency differences >20% compared to East Asian-specific allele frequencies. The genotype scaffold for each study was then imputed to the 1000G Phase 1 or 3 reference panel[53] using minimac3[54] or IMPUTEv2[55]. In BMI-unadjusted analyses, all studies were imputed to 1000G Phase 3. In BMI-adjusted and sex-stratified analyses, all studies were imputed to 1000G Phase 3 except for a subset of Biobank Japan[17], which was imputed to the 1000G Phase 1 reference panel.

Study-level association analyses

Within each study, all variants were tested for association with T2D assuming an additive model of inheritance within a regression framework, including age, sex, and other study-specific covariates (Supplementary Table 3). To account for population structure and relatedness, association analyses were either performed using FIRTH[56] or mach2dat with additional adjustment for principal components in unrelated individuals or a linear mixed model with kinship matrix implemented in BOLT-LMM[57]. In studies analyzed with the linear mixed model, allelic effects and standard errors were converted to the log-odds scale that accounts for case-control imbalance[58]. Within each study, variants were removed if i) imputation quality score was poor (minimac3 r2<0.30; IMPUTE2 info score <0.40); ii) combined case control minor allele count <5; or iii) standard error of the log-OR>10. For a subset of the studies, BMI was added as an additional covariate, and association analyses were also performed separately in males and females. For each study and model, association statistics were corrected with genomic control inflation factor[59] calculated from common variants (MAF≥5%) (Supplementary Table 3). For BBJ, we applied the genomic control inflation factor 1.21 as reported[2].

Sex-combined meta-analysis

We combined study-level association statistics using fixed-effects meta-analysis with inverse-variance weighting of log-ORs implemented in METAL[60]. Variants with allele frequency differences >20% between 1000G Phase 1 and 3 panels were excluded from the meta-analysis. To assess excess inflation arising from cryptic relatedness and population structure, we applied LD score regression to the meta-analysis summary statistics to estimate residual inflation of summary statistics, using a set of 1,889 unrelated Chinese individuals from the Singapore Chinese Eye Study[61]. The LD score regression intercepts were 0.993 for BMI-unadjusted, and 1.0163 for BMI-adjusted models. As the LD score regression intercepts indicated absence of excess inflation, the meta-analysis results were corrected for inflation using these LD score regression intercepts. For subsequent analyses, we considered only variants that were present in at least 50% of the effective sample size Neff [computed as 4/(1/Ncases + 1/Ncontrols)][60]. Heterogeneity in allelic effect sizes between studies were assessed with fixed-effects inverse-variance weighted meta-analysis P. We further compared the genetic effects from BMI-unadjusted and BMI-adjusted models using fixed-effects inverse-variance weighted meta-analysis P. Loci were defined as novel if the lead variant is: (1) at least 500 kb away and confirmed by GCTA to be conditionally independent from previously reported T2D-associated variants in any ancestry, and (2) assessed using LocusZoom plots and detailed literature review to be away from known loci with extended LD. Lead variants mapping to loci already associated with other glycemic traits were still considered novel for the association with T2D.

BMI adjustment analyses and effect size comparison

For the subset of studies with both BMI-unadjusted and BMI-adjusted models, we compared the effect sizes and heterogeneity of the lead variants using the standardized mean difference to account for the correlation between the two models[1]. We calculated the Pearson correlation coefficient between effect sizes from the BMI-unadjusted and BMI-adjusted models for all 13.2M variants genome-wide (r=0.98) and for the lead variants at 189 T2D-associated loci (r=0.99).

Sex-differentiated meta-analysis

The meta-analyses described above were repeated for males and females separately. The male-specific meta-analyses included up to 28,027 cases and 89,312 controls (Neff = 65,660) and the female-specific analyses included up to 27,370 cases and 135,055 controls (Neff = 70,332). LD score regression intercepts were 1.0044 for BMI-unadjusted and 1.0045 for BMI-adjusted models in males and 1.0050 for BMI-unadjusted and 1.0187 for BMI-adjusted models in females. We further performed a test for heterogeneity in allelic effects between males and females as implemented in GWAMA[62,63].

Detection of distinct association signals

To detect multiple distinct association signals at each associated locus, we combined overlapping loci when the distance between any pair of lead variants was <1 Mb. We then performed approximate conditional analyses using GCTA[12] with genome-wide meta-analysis summary statistics and LD estimated from 78,000 samples from the Korean Biobank Array[52]. We note the limitations in using a single population reference panel for LD estimation for a meta-analysis of diverse East Asian populations. We present all distinct signals at conditional threshold of P<1x10−5, but we suggest that readers exhibit caution and limit inferences from these analyses to signals that show the strongest evidence of association.

Comparing loci effects between East Asian and European populations

We compared the per-allele effect sizes of lead variants identified from the East Asian BMI-unadjusted sex-combined meta-analysis (183 lead variants) and European BMI-unadjusted sex-combined meta-analysis[1] (231 lead variants; Supplementary Table 7). Across the 414 associated variants from the two ancestries, 12 lead variants overlapped, resulting in 402 unique variants. As the variants in the European analysis were imputed using the Haplotype Reference Consortium reference panel and did not include indel variants, a variant in strong LD (East Asian r2>0.90) with the lead East Asian variant was used when the lead variant was an indel, when possible. If the lead East Asian variant or a variant in strong LD (East Asian r2>0.90) was not available in the European data from DIAMANTE, we used results from a previous European type 2 diabetes meta-analysis[64]. The effect size comparison plot was restricted to 332 variants where data was available for both ancestries (Figure 2A). For loci that were significant in both the East Asian and European meta-analyses, if the lead variants were different, both lead variants were plotted (see Supplementary Table 7). Effect size plots were further restricted to: i) 278 lead variants with MAF≥5% in both East Asian and European meta-analyses (Extended Data Figure 7); ii) 203 lead variants significant in the East Asian meta-analysis (Extended Data Figure 8A); and iii) 234 lead variants significant in the European meta-analysis (Extended Data Figure 8B). Differences in effect sizes between the two populations could be due to differences in imputation quality with different reference panels.

Associations with other metabolic traits and outcomes

We examined publicly available GWAS summary statistics (mostly available through the Type 2 Diabetes Knowledge Portal[37]) to explore associations of the lead variant at the 92 loci for which there are no genome-wide significant European variants within 500 kb (listed in Supplementary Table 9). Association statistics from the following consortia were available for query on the portal (last accessed August 28, 2019): coronary artery disease from CARDIoGRAM[65], BMI and waist-hip-ratio from GIANT[66,67], lipid traits from GLGC[68], and glycemic traits from MAGIC[69-73]. Additionally, we used available data from AGEN East Asian meta-analyses for lipids[74], along with the phenotypic data from the UK Biobank[75], BioBank Japan[24,38], and blood pressure data from ICBP[76]. For this analysis, we looked up the effect size and P-value of the East Asian lead variants in the other datasets. If the variant-trait association reached at least nominal significance (P<1x10−3), we included the lookup results in Supplementary Table 10. When the lead East Asian variant was missing in the prior GWAS data, we reported it as “NF” (not found) in the table.

Colocalization with expression quantitative trait loci (eQTL)

We searched publicly available eQTL databases such as GTEx v7[77] and the Parker lab Islet Browser[19], to identify cis-eQTLs at the novel loci in adipose (subcutaneous and visceral), blood, pancreas, pancreatic islet, and skeletal muscle tissue. We also searched for cis-eQTLs in subcutaneous adipose tissue data from the METSIM study[6], whole blood[78], and peripheral blood (BioBank Japan; http://jenger.riken.jp/en/result). Colocalized eQTLs were identified if the lead expression level-associated variant and the GWAS lead variant were in high LD (r2>0.80) in Europeans to accommodate the predominantly European eQTL data. Reciprocal conditional analyses were also performed using the METSIM data to determine if the GWAS lead variant and the lead expression-associated variant were part of the same eQTL signal.

Literature review

We conducted a traditional literature review to identify candidate genes at each novel locus using NCBI Entrez Gene, PubMed and OMIM. We included gene symbols and the following keywords as search terms in PubMed: diabetes, glucose, insulin, islet, adipose, muscle, liver, obesity. A gene was considered a potential candidate if an apparent link to T2D biology existed based on prior studies of gene function.

Functional annotation and experimentation at NKX6-1

We used ENCODE[79], ChromHMM[80], and Human Epigenome Atlas[81] data available through the UCSC Genome Browser to identify candidate variants at the association signal near NKX6-1 that overlapped open-chromatin peaks, ChromHMM chromatin states, and chromatin-immunoprecipitation sequencing (ChIP-seq) peaks of histone modifications H4K4me1, H3K4me3, and H3K27ac, and transcription factors in the pancreas and pancreatic islets. MIN6 mouse insulinoma cells (FROM ATCC)[82] were cultured in DMEM (Sigma) supplemented with 10% FBS, 1mM sodium pyruvate, and 0.1 mM beta-mercaptoethanol. The cell cultures were maintained at 37° C with 5% CO2. To measure variant allelic differences in enhancer activity at the NKX6-1 locus, we designed oligonucleotide primers (forward: CCCTAGTAATGCCCTTTTTGTT; reverse: TCAGCCTGAGAAGTCTGTGA) with KpnI and Xhol restriction sites, and amplified the 400-bp DNA region (GRCh37/hg19 -chr4: 85,339,430-85,339,829) around rs117624659. As previously described[80], we ligated amplified DNA from individuals homozygous for each allele into the multiple cloning site of the pGL4.23 (Promega) minimal promoter luciferase reporter vector in both the forward and reverse orientations with respect to the genome. Clones were isolated and sequenced for genotype and fidelity. 2.1x105 MIN6 cells were seeded per well and grown to 90% confluence in 24-well plates. We co-transfected five independent luciferase constructs and Renilla control reporter vector (phRL-TK, Promega) using Lipofectamine 2000 (Life Technologies) and incubated. 48-hours post-transfection, the cells were lysed with Passive Lysis Buffer (Promega). Luciferase activity was measured using the Dual-luciferase Reporter Assay System (Promega) per manufacturer instructions and as previously described[83]. MIN6 cell lines were authenticated through genotyping and tested negative for mycoplasma contamination.