A full description of the methods used in this study is available as Supplementary Methods.

Data reporting

The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Sample selection

We drew samples for exome sequencing from six consortia, most of which consisted of multiple studies and are described fully in Supplementary Table 1. T2D case status was determined according to study-specific criteria described in full in in Supplementary Table 1 and the Supplementary Methods. All individuals provided informed consent and all samples were approved for use by their institution’s institutional review board or ethics committee, as previously reported10,46,47,48. Samples that were newly sequenced at The Broad Institute as part of T2D-GENES, SIGMA and ProDiGY are covered under Partners Human Research Committee protocol 2017P000445/PHS ‘Diabetes Genetics and Related Traits’.

Data generation

The details of data generation, variant calling, quality control and variant annotation are described in full in the Supplementary Methods. In brief, for each consortium, sequencing data were aggregated (if previously available) or newly generated (if not) and then processed through a standard variant calling pipeline. We then measured samples and variants according to several metrics indicative of sequencing quality, excluding those that were outliers relative to the global distribution (Supplementary Fig. 1, Supplementary Table 2). These exclusions produced a ‘clean’ dataset of 49,484 samples and 7.02 million variants.

Following initial sample and variant quality control, we performed additional rounds of sample exclusion from association analysis (Extended Data Fig. 2). We also excluded the 3,510 childhood diabetes cases from the SEARCH and TODAY studies based on an analysis that suggested their lack of matched controls would induce artefacts in gene-level association analyses (Supplementary Fig. 17). These exclusions produced an ‘analysis’ dataset of 45,231 individuals and 6.33 million variants. A power analysis of this dataset is presented in the Supplementary Methods.

After these three rounds of sample exclusions, we estimated—within each ancestry—pairwise identity-by-descent values, genetic relatedness matrices and principal components for use in downstream association analyses. We used the identity-by-descent values to generate lists of unrelated individuals within each ancestry, excluding 2,157 individuals to produce an ‘unrelated analysis’ set of 43,090 individuals (19,828 cases and 23,262 controls) and 6.29 million non-monomorphic variants. We used this set of individuals and variants for single-variant and gene-level tests (described below) that required an unrelated set of individuals.

We annotated variants with the ENSEMBL Variant Effect Predictor49 (VEP, version 87). We produced both transcript-level annotations for each variant as well as a ‘best guess’ gene-level annotation using the –flag-pick-allele option (with ranked criteria described in the Supplementary Methods). We used the VEP LofTee (https://github.com/konradjk/loftee) and dbNSFP (version 3.2)50 plugins to generate additional bioinformatics predictions of variant deleteriousness; from the dbNSFP plugin, we took annotations from 15 different bioinformatics algorithms (listed in Extended Data Fig. 5) and then added annotations from the mCAP51 algorithm. As these annotations were not transcript-specific, we assigned them to all transcripts for the purpose of downstream analysis.

Although we incorporated both transcript-level and gene-level annotations into gene-level analyses (see below), all single-variant analyses reported in the manuscript or figures are annotated using the ‘best guess’ annotation for each variant.

Single-variant association analysis in sequencing data

To perform single-variant association analyses, we first stratified samples by cohort of origin and sequencing technology (with some exceptions described in the Supplementary Methods), yielding 25 distinct sample subgroups (Extended Data Fig. 3). For each subgroup, we performed additional variant quality control beyond that used for the ‘clean’ dataset, excluding variants according to subgroup-specific criteria described in Extended Data Fig. 3; in general, these criteria were strict—particularly for multiallelic variants and X-chromosome variants. We verified that these filters led to a well-calibrated final analysis through inspection of quantile–quantile plots within and across ancestries (Extended Data Fig. 4).

For each of the 25 sample subgroups, we then conducted two single-variant association analyses: one of all (including related) samples using the (two-sided) EMMAX test52 and one of unrelated samples using the (two-sided) Firth logistic regression test53. Both analyses included covariates for sequencing technology, and the Firth analysis included covariates for principal components of genetic ancestry (those among the first 10 that showed P < 0.05 association with T2D).

We then conducted a 25-group fixed-effect inverse-variance weighted meta-analysis for each of the Firth and EMMAX tests, using METAL54. We used EMMAX results for association P values and Firth results for effect size estimates.

Additional analysis of rs145181683

To assess whether the rs145181683 variant in SFI1 (P = 3.2 × 10−8 in the exome-sequencing analysis) represented a true novel association, we obtained association statistics from 4,522 Latinos55) who did not overlap with the current study. On the basis of the odds ratio (1.19) estimated in our analysis and the MAF (12.7%) in the replication sample, the power was 91% to achieve P < 0.05 under a one-sided association test. The observed evidence (P = 0.90, odds ratio = 1.00) did not support rs145181683 as a true T2D association. Further investigation of this lack of replication evidence suggested that, although the association from our sequence analysis is unlikely to be a technical artefact (genotyping quality was high), it could possibly be a proxy for a different (Native American-specific) non-coding causal variant (full details are available in the Supplementary Methods). Further fine-mapping and replication efforts will be necessary to test this hypothesis.

Gene-level analysis

For each gene, following previous studies10,56,57, we separately tested seven different ‘masks’ of variants grouped by similar predicted severity (defined in Extended Data Fig. 5). For each gene and each mask, we created up to three groupings of alleles, corresponding to different transcript sets of the gene; for many genes, two or more of these allele groupings were identical.

Before running gene-level tests, we performed additional quality control on sample genotypes. For each of the 25 sample subgroups (the same as used for single-variant analyses), we identified variants that failed subgroup-specific quality control criteria (shown in Extended Data Fig. 5) and set genotypes for these variants in all individuals in the subgroup as ‘missing’.

We conducted two gene-level association tests: a burden test, which assumes all analysed variants within a gene are of the same effect, and SKAT15, which allows variability in variant effect size (and direction); each of these tests is two-sided. We performed each test across all unrelated individuals with 10 principal components of genetic ancestry, sample subgroup and sequencing technology as covariates. As this ‘mega-analysis’ strategy was different from the meta-analysis strategy that we used for single-variant analyses, as a quality control exercise we conducted a single-variant mega-analysis and found that its results showed broad correlation with those from the original meta-analysis (Supplementary Fig. 18).

We then developed two methods to consolidate the 2 × 7 = 14 P values produced for each gene (described in full in Extended Data Fig. 5, Supplementary Methods and Supplementary Figs. 5, 6). First, we corrected the smallest P value for each gene according to the effective number of independent masks tested for the gene (variable, but on average 3.6), based on the gene-specific correlation of variants across masks58 (referred to as the minimum P-value test; Supplementary Fig. 19). Second, we tested all nonsynonymous variants (that is, missense, splice site and protein-truncating mutations), but weighted each variant according to its estimated probability of causing gene inactivation9 (referred to as the weighted test, which essentially assessed the effect of gene haploinsufficiency from combined analysis of protein-truncating and missense variants; Supplementary Fig. 6). We verified that these two consolidation methods were well-calibrated (Extended Data Fig. 6) and broadly consistent yet distinct: across the 10 most significantly associated genes, P values were nominally significant using both methods for 8 genes but varied by 1–3 orders of magnitude (Extended Data Table 2).

Because each gene mask could in fact represent up to three sets of alleles (owing to the transcript-specific annotation strategy that we used), for each of the four analyses multiple P values were possible for some genes. To produce a single gene-level P value for each of the four analyses, we thus collapsed (for each gene) the set of P values across transcript sets into a single gene-level P value using the minimum P-value test.

We used a conservative Bonferroni-corrected gene-level exome-wide significance threshold of P = 0.05/(2 tests × 2 consolidation methods × 19,020 genes) = 6.57 × 10−7. For each gene referenced in the manuscript, we report the P value and odds ratio from the analysis that achieved the lowest P value for the gene.

Gene-level analysis near T2D GWAS signals

In principle, a nearby common-variant association could lead to over- or underestimation of the strength of a gene-level association59. To assess whether differential patterns of rare variation across common-variant haplotypes could significantly affect our gene-level results, we conducted two analyses (described in the Supplementary Methods) and found no evidence that confounding from common-variant haplotypes was primarily responsible for the associations that were observed in our gene-level analyses.

Further exploration of significant gene-level associations

For our exome-wide significant gene-level associations (MC4R, PAM and SLC30A8), we conducted additional gene-level analyses to dissect the aggregate signals that were observed. First, we performed tests by progressively removing alleles in order of lowest single-variant analysis P value, in order to understand the (minimum) number of alleles that contributed statistically to the aggregate signal. Second, we performed tests conditional on each allele in the sequence (that is, calculating separate models with each individual allele as a covariate), and we then compared the resulting P values to the full gene-level P value, in order to assess the contribution of each allele individually to the signal. Finally, for MC4R, we conducted an analysis with an added sample covariate for body-mass index and found that it, as shown previously60,61, reduces the significance of both the Ile269Asn single-variant signal (P = 1.0 × 10−5) and the gene-level signal not attributable to Ile269Asn (P = 0.035).

To evaluate which ancestries contributed variants to MC4R, SLC30A8, and PAM, we calculated the proportion of variants in each signal unique to an ancestry and also compared the significance and direction of effect of each signal across ancestries. Across the three signals, 68.4% (287 out of 419) of variants in total were unique to one ancestry (63.9% for MC4R, 67.0% for SLC30A8 and 71.6% for PAM). Each signal had a direction of effect that was consistent across all five ancestries and each signal achieved P < 0.05 in at least two ancestries (MC4R in East-Asians and Hispanics; SLC30A8 in all ancestries other than African-Americans; and PAM in Europeans, South-Asians and Hispanics).

Analysis of exomes from the Geisinger Health System

We obtained gene-level association results that were previously computed from an analysis of 49,199 individuals (12,973 T2D cases and 36,226 controls) from the Geisinger Health System (GHS). Association statistics were available for 44 out of the 50 genes with the strongest gene-level associations from our study. A power analysis of the GHS replication analysis is available in the Supplementary Methods.

GHS sequencing data were processed and analysed as previously described24, and variants were grouped into four (nested) masks (roughly corresponding to the LofTee, 5/5, 1/5 1% and 0/5 1% masks; more details are available in the Supplementary Methods). For each mask, association results were computed using two-sided logistic regression under an additive burden model (with phenotype regressed on the number of variants carried by each individual) with age, age2 and sex as covariates. To produce a single GHS P value for each gene, we applied the minimum P-value procedure across the four mask-level results.

Analysis of exomes from the CHARGE consortium

We collaborated with the CHARGE consortium to analyse the 50 genes with the strongest gene-level associations from our study in 12,467 individuals (3,062 T2D cases and 9,405 controls) from their previously described study62,63. A power analysis of the CHARGE replication analysis is available in the Supplementary Methods.

Variants in the CHARGE exomes were annotated and grouped into seven masks using the same procedure as for the original exome-sequencing analysis. Burden and SKAT association tests were then performed in the Analysis Commons64 using a two-sided logistic mixed model65 assuming an additive genetic model and adjusted for age, sex, study, race and kinship. To produce a single CHARGE P value for each gene, we applied the minimum P-value procedure across the seven mask-level results, as for the GHS analysis.

Meta-analysis with CHARGE and GHS

We conducted a meta-analysis among our original burden analysis and those of CHARGE and GHS. For each gene, we selected the mask that achieved the lowest P value in our original analysis and conducted a two-sided sample-size weighted meta-analysis with the results from CHARGE and GHS for the same mask (or an analogous mask as defined in the Supplementary Methods).

Investigation of the UBE2NL association

We investigated the novel association that was found in the gene-level meta-analysis (UBE2NL, meta-analysis P = 5.6 × 10−7) in more detail. The UBE2NL burden signal was due to five PTVs in the original analysis (observed in 29 cases and 1 control; all of which had high (>45×) sequencing coverage; Supplementary Table 8) and was replicated at P = 0.02 in CHARGE; UBE2NL results were not available in GHS. As UBE2NL lies on the X chromosome, we conducted a sex-stratified analysis of the original samples and observed independent associations in both men (P = 5.7 × 10−4) and women (P = 1.6 × 10−3). UBE2NL does not lie near any known GWAS associations (http://www.type2diabetesgenetics.org/) and has few available references66,67,68, suggesting that it may be a novel T2D-relevant gene, although further replication will be important to establish its association.

Evaluation of directional consistency between exome-sequencing, CHARGE and GHS analyses

We examined the concordance of direction of effect size estimates (that is, both odds ratios of >1 or <1) between burden tests from our original exome-sequencing analysis and those from CHARGE and GHS. For the 46 genes advanced for replication with burden P < 0.05 for at least one mask (that is, ignoring those with evidence for association only under the SKAT model), we compared the direction of effect estimated for the mask with lowest P-value mask to that estimated for the same (or analogous) mask in the GHS or CHARGE analysis. We then conducted a one-sided exact binomial test to assess whether the fraction of results with consistent direction of effects was significantly greater than expected by chance.

Gene set analysis in sequencing data

We curated 16 sets of candidate T2D-relevant genes, defined in Supplementary Table 9 with criteria as specified in the Supplementary Methods. For each gene set, we constructed sets of matched genes with similar numbers and frequencies of variants within them (details are provided in the Supplementary Methods). A sensitivity analysis of this matching strategy is presented in the Supplementary Methods.

To conduct a gene set analysis, we then combined the genes in the gene set with the matched genes. Within the combined list of genes, we ranked genes using the P values observed for the minimum P-value burden test. We then used a one-side Wilcoxon rank-sum test to assess whether genes in the gene set had significantly higher ranks than the comparison genes.

Use of gene-level associations to predict effector genes

To assess whether gene-level associations from exome sequencing—which are composed mostly of rare variants independent of any GWAS associations—could prioritize potential effector genes within known T2D GWAS loci, we first assessed whether predicted effector genes (based on common-variant associations) were also enriched for rare coding variant associations. Our analysis (described in full in the Supplementary Methods) indicated that effector genes predicted from common coding variant associations do show significant enrichments (P = 8.8 × 10−3), but effector genes predicted from transcript-level associations do not (P = 0.72).

We then curated a list of 94 T2D GWAS loci, and 595 genes that were within 250 kb of any T2D GWAS index variant, from a 2016 T2D genetics review69 and observed 40 with a P < 0.05 gene-level signal (Supplementary Table 12), greater than the 595 × 0.05 = 29.75 expected by chance (P = 0.038). Only three (SLC30A8, PAM and HNF1A) were from the list that we curated of 11 genes with causal common coding variants6. We found that these 40 genes were significantly more enriched for protein interactions (P = 0.03; observed mean = 11.4, expected mean = 4.5) than the 184 genes implicated based on proximity to the index SNP (P = 0.64; observed mean = 21.1, expected mean = 21.9), although evaluation of the biological candidacy of these genes will ultimately require in-depth functional studies70. Rare coding variants could therefore, in principle, complement common-variant fine-mapping71,72 and experimental data4,70 to help to interpret T2D GWAS associations; however, our results indicate that much larger sample sizes and/or orthogonal experimental data will be required to clearly implicate specific effector genes. A full description of this analysis is included in the Supplementary Methods.

Use of gene-level associations to predict direction of effect

To assess whether gene-level association analyses of predicted deleterious variants could be used to predict therapeutic direction of effect, we compared odds ratios estimated from a modified weighted burden test procedure (described in the Supplementary Methods) to those expected for T2D drug targets (assuming agonist targets to have true odds ratios > 1 and inhibitors to have true odds ratios < 1). For a similar comparison to expectations for mouse gene knockouts, we used the relationship between mouse phenotype and human phenotype specified in the Supplementary Methods. Genes present in two gene sets with opposite expected direction of effects were excluded from this analysis.

Collection and analysis of SNP array data

To compare discoveries from our exome-sequencing analyses to discoveries possible from common-variant GWAS of the same samples, we aggregated all available SNP array data for the exome-sequenced samples (18,233 cases and 17,679 controls; Supplementary Table 13). After sample and variant quality control (described in the Supplementary Methods), we imputed variants from the 1000 Genomes Phase 332 (1000G) and Haplotype Reference Consortium33 (HRC) reference panels using the Michigan Imputation Server73. We used 1000G-based imputation for all association analyses and HRC-based imputation to assess the number of exome-sequence variants imputable from the largest available European reference panel (details available in the Supplementary Methods).

After imputation, we performed sample and variant quality control, as well as two-sided association tests, analogous to the exome-sequence single-variant analyses. In contrast to the exome-sequencing analyses, a quantile–quantile plot suggested that the associations from the EMMAX test were not well calibrated, and we therefore used only the Firth test (that is, for both P values and odds ratios) in the imputed GWAS analysis.

To conduct gene set analysis with the imputed GWAS data, we first used the method implemented in MAGENTA74 to calculate gene scores from the imputed GWAS single-variant association results. Following the same protocol as for gene set analysis from the exome-sequencing results, we then conducted a one-sided Wilcoxon rank-sum test to compare the gene scores to those of matched comparison genes. We followed the same approach for the gene set analysis that we conducted in a larger, previously published13 GWAS.

LVE calculations

To calculate LVEs, we used a previously presented formula75 (equations are available in the Supplementary Methods) to calculate the LVE of a variant with three genotypes (AA, Aa and aa) and corresponding relative risks (1, RR 1 and RR 2 ). When presenting the strongest LVE values for the imputed GWAS analysis, we only considered variants that were genotyped in at least 10,000 individuals to avoid potential artefacts that result from a spurious association in a small-sample subgroup. For gene-level LVE calculations, we used the variant mask with lowest P value to calculate LVEs. We also conducted a sensitivity analysis to bound the extent to which our gene-level LVE estimates might be biased downwards due to their inclusion of benign alleles; this analysis (described in full in the Supplementary Methods) produced upper bounds of gene-level LVEs that were at most twofold higher than the point estimates.

Prediction of LVE explained by the top 100 and top 1,000 gene-level associations

To forecast the LVE that will be explained once 100 (or 1,000) significant T2D gene-level associations are detected, we applied a previously suggested model34 in which the LVE of a gene is related to its rank in the overall gene-level P-value distribution. Specifically, the model is LVE n = ean + b where LVE n is the LVE of the gene with nth lowest gene-level P value. We fitted this model using linear regression to the top 50 genes in our analysis (Supplementary Fig. 20), yielding estimates of a = −0.044 and b = −7.07. We then calculated the LVE of the top 100 (or 1,000) genes by summing the actual LVE of the top three signals (which achieved exome-wide significance in our analysis) with the LVE predicted by the model for genes ranked 4–100 (or 4–1,000).

Estimated power to detect gene-level associations with T2D drug targets

To estimate the power of future studies to detect gene-level associations in genes with effect sizes similar to those for established T2D drug targets, we used aggregate allele frequencies and odds ratios estimated from our gene-level analysis and an assumed prevalence of K = 0.08 to calculate a proxy for true population frequencies and relative risks. For each gene, we used odds ratios and frequencies from the variant mask that yielded the strongest gene-level association. Because, on average, these drug targets had five effective tests per mask, we used an exome-wide significance threshold of α = 1.25 × 10−7 for power calculations. We calculated power as previously described76.

The ranges given in the main text (75,000–185,000 disease cases) represent the numbers from the power calculations for INSR (the drug target with the highest observed effect size) and IGF1R (the drug target with the lowest observed effect size other than KCNJ11 and ABCC8). We excluded KCNJ11 and ABCC8 from this reported range, given that a mixture of risk-increasing and risk-decreasing variants in these genes probably diluted their burden signals. We did not account for uncertainty in estimated odds ratios or aggregate variant frequency in these calculations, as no genes had 95% confidence intervals that that did not overlap odds ratio = 1.

Interpretation of suggestive associations

We quantified the PPA for nonsynonymous variants observed in our dataset as a function of association strength measured by single-variant P values. We define a true association as a variant that, when studied in larger sample sizes, will eventually achieve statistical significance owing to a true odds ratio ≠ 1. We distinguish true associations from causal associations: causally associated variants are the subset of truly associated variants in which the variant itself is causal for the increase in disease risk, as opposed to being truly associated due to linkage disequilibrium (LD) with a different causally associated variant (that is, an ‘LD proxy’). An overview of the method that we developed for PPA calculations is provided in Extended Data Fig. 7, and a full description of the method is included in the Supplementary Methods. Here, we outline the steps in the approach.

First, for various single-variant P-value thresholds in the exome-sequencing analysis, we calculated the fraction of variants that reached this threshold with directions of effect concordant with those of an independent exome array study10. For example, 61.3% of nonsynonymous variants within T2D GWAS loci that reached P < 0.05 in the exome-sequencing analysis had concordant directions of effect with the independent study, a fraction that decreased (as expected) for higher P-value thresholds (for example, 49.4% at P > 0.5) or when only variants outside of T2D GWAS loci were analysed (51.9% at P < 0.05).

Second, we derived an equation to convert the fraction of concordant associations to an estimated proportion of true associations. This value provides a PPA estimate, as a function of P value, for an arbitrary variant in the set initially used to calculate direction of effect concordances. We computed separate mappings for arbitrary nonsynonymous variants (using all exome-wide nonsynonymous variants) and one for nonsynonymous variants within GWAS loci (using only nonsynonymous variants within the 94 T2D GWAS loci). We note that the mapping produced from our analysis applies only to the results from the current study: because other studies have different sample sizes and may apply different statistical tests, the mapping would need to be recomputed to interpret the associations of other studies using the same method.

Third, we converted PPA estimates to estimates of the posterior probability of causal associations (PPA c ). This conversion requires estimates of the fraction of coding variant associations that are causal (as opposed to LD proxies). We explored several values for this parameter, as described in the Supplementary Methods and shown in Extended Data Fig. 8.

Fourth, we extended PPA estimates to incorporate gene-specific priors by mapping posterior odds of causal association (PO c ) to a Bayes factor for causal association (BF c ). This calculation requires a set of training variants with a known prior. For this training set, we use nonsynonymous variants within GWAS loci and modelling assumptions for their prior. Details of this model are described in the Supplementary Methods and a sensitivity analysis of its assumptions is shown in Extended Data Fig. 8.

Finally, as a preliminary estimate of a principled prior likelihood for genes in the mouse NIDD gene set, we estimated the proportion of non-null associations across all genes in the set. To use true prior data (rather than associations from the current study), we calculated gene-level P values for each gene in the set using the MAGENTA74 algorithm applied to a recent transethnic T2D GWAS13. We then used a previously developed approach40,77 that models the distribution of observed P values as a mixture of uniform (representing the null distribution) and beta (representing the non-null distribution) distributions, yielding a prior value of 23.2%.

Our PPA c calculations currently have several limitations. They apply only to single-variant associations and not (yet) to gene-level associations; extending them to apply to gene-level associations would avoid the possibility of conflicting results among variants within a gene but would require larger-scale gene-level replication data than that we had available in the current analysis. Additional work will also be needed to generate data and develop methods to estimate objective rather than subjective gene priors (researchers can often overestimate evidence of disease relevance for genes in which they have invested considerable effort), to reduce dependence of our conclusions on modelling assumptions (Extended Data Fig. 8) and to explore the extent to which the large number of variant associations that we predict from our data localize to specific gene or variant functional annotations78.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.