Study participants and ethics approval
Our MWAS for incident T2D involves the use of data from ten prospective cohorts, including the Nurses’ Health Study (NHS; initiated in 1976 with 121,701 female nurses aged 30–55 years9,50), NHS2 (started in 1989 with 116,429 female nurses aged 25–42 years9,50), Health Professionals Follow-Up Study (HPFS; started in 1986 with 51,529 male health professions aged 40–75 years9), Hispanic Community Health Study/Study of Latinos (SOL; enrolled 16,415 Hispanic/Latino adults aged 18–74 years during 2008–201151,52), Women’s Health Initiative (WHI; initiated in 1993 enrolling 68,132 women aged 50–79 years to one of three clinical trials or an observational study53), Atherosclerosis Risk in Communities (ARIC) study (enrolled 15,792 mostly Black and white US adults aged 45–64 years during 1987–198954), Framingham Heart Study Offspring cohort (FHS; enrolled 5,124 adults; we focused on those attended the fifth examination during 1991–1995), Multi-Ethnic Study of Atherosclerosis (MESA; initiated in 2000 with 6,814 adults aged 45–84 years55,56), the Boston Puerto Rican Health Study (BPRHS; enrolled 1,500 self-identified Puerto Rican adults aged 45–75 years) and the Prevención con Dieta Mediterránea Study (PREDIMED; a 5-year dietary trial with 7,447 adults aged 55–80 years57). In each cohort, comprehensive data on demographics, medical and family history, diet, lifestyle and other health information were collected at baseline and were updated during longitudinal follow-ups. Blood samples were collected at baseline and/or during follow-ups. Our MWAS for incident T2D included participants with qualified metabolomics data, and were free of diabetes, cardiovascular disease and cancer at study baseline. The final analysis included 6,890 participants from NHS; 3,692 from NHS2 and 2,529 from HPFS; 2,821 from SOL; 1,392 from WHI; 1,288 white and 1,433 Black participants from ARIC; 1,424 from FHS; 902 from MESA; 378 from BPRHS and 885 from PREDIMED (Extended Data Table 1). Each study was approved by Institutional Review Boards at respective institutions or study centers, and all participants provided informed consent. Our GWAS for metabolites included participants from eight cohorts comprising NHS, NHS2, HPFS, SOL, WHI, ARIC, FHS and, in addition, the Cardiovascular Health Study (CHS; enrolled 5,201 adults during 1989–1990 and 678 predominantly Black participants in 1992–199358,59) (Supplementary Table 7). The detailed descriptions of the design, data collection, ethical review of each cohort, and our inclusion and exclusion criteria are provided in Supplementary Methods.
Ascertainment of T2D
In all cohorts, incident T2D was defined when a participant was free of diabetes at baseline but was identified as having T2D during longitudinal follow-up. Detailed information on diagnosis criteria in each cohort is included in Supplementary Methods, and follow-up years and numbers of incident cases are listed in Extended Data Table 1. Briefly, in NHS/HPFS, T2D were identified by follow-up questionnaires, and confirmed through a supplementary questionnaire based on diagnostic criteria from the National Diabetes Data Group before 199860 and the American Diabetes Association (ADA) criteria after 199861,62. In SOL, T2D was defined if a participant had fasting glucose ≥7.0 mmol l−1, fasting ≤8 h and nonfasting glucose ≥11.1 mmol l−1, post oral glucose tolerance test glucose ≥11.1 mmol l−1, HbA1c ≥ 6.5%, current use of antidiabetic medications or self-reported physician-diagnosed diabetes63. In WHI, T2D was determined based on self-reported history of diabetes or using antidiabetic medications (pills or shots) in any visits/interviews. In ARIC and FHS, T2D was diagnosed if a person had fasting glucose ≥7.0 mmol l−1, fasting ≤8 h and nonfasting glucose ≥11.1 mmol l−1, or current use of antidiabetic medications with ARIC further considering self-reported physician-diagnosed diabetes64,65. T2D cases in MESA and BPRHS were determined according to the ADA criteria66, which included fasting plasma glucose level ≥7.0 mmol l−1 or the use of antidiabetic medications or insulin56,67. In PREDIMED, T2D was adjudicated through blind assessment by a Clinical Endpoint and Adjudication of Events Committee, based on the ADA criteria68.
Assessment of diet, lifestyle factors and covariates
Detailed information on data collection in each cohort is in Supplementary Methods. Briefly, demographic factors (for example, self-reported sex, and race and ethnicity), socioeconomic status, health information (for example, medical conditions and family history) and lifestyle (for example, smoking history and PAs), anthropometrics and blood pressure, were collected at baseline and follow-up visits, through self-administrated questionnaires, or in-person or telephone-based interviews by trained staff. PA was quantified as metabolic equivalent (MET) in hours per week. We calculated BMI based on baseline weight and height, and WHR based on waist and hip circumferences. Blood clinical biomarkers were measured using standard assays. Among participants with serum creatinine data, eGFR was estimated using the CKD-EPI (Chronic Kidney Disease Epidemiology Collaboration) formula, based on age, sex and race in NHS/HPFS, WHI, ARIC and PREDIMED69, and standard reference equations for Hispanics adjusting for age and sex in SOL. In PREDIMED, two propensity scores were estimated to account for the probability of assignment to intervention groups57.
In NHS/HPFS, diet was assessed using a semi-quantitative food frequency questionnaire (FFQ) every 4 years; in our analysis we averaged the intakes from the two FFQs closest to the time of blood draw (NHS: 1986 and 1990; NHS2: 1995 and 1999; HPFS: 1994 and 1998). In WHI, ARIC, FHS, MESA and BPRHS, diet was similarly assessed by FFQs designed and validated for application to their targeted populations (for example, multiethnic and geographically diverse populations in WHI70,71,72 and Puerto Rican population in BPRHS73). In SOL, diet was assessed using two 24-h dietary recalls and a food propensity questionnaire74. The overall dietary quality was assessed by the Alternate Healthy Eating Index-2010 (AHEI-2010)75 in all cohorts except for the PREDIMED trial, in which it was assessed by a 14-item Mediterranean Diet Adherence Screener score57. In NHS/HPFS, SOL and WHI, we also calculated baseline consumptions of 15 main food groups in the unit of servings per day.
Metabolomic profiling, quality control and data harmonization
Metabolomic profiling in NHS/HPFS, WHI, MESA, PREDIMED, FHS and CHS was conducted with the Metabolomics Platforms at the Broad Institute of MIT and Harvard University, using three to four complementary LC–MS methods9,65,76. Metabolomic profiling in SOL and ARIC (serum samples) and BPRHS (plasma samples) was conducted using LC–MS based methods by the Metabolon DiscoveryHD4 Panel at the Metabolon Inc.63,77,78. Detailed protocols for both platforms have been described previously53,79.
Data processing was conducted within each study and, if applicable, separately within each batch (or substudy) if several batches/substudies were conducted within a cohort. Samples were removed if their metabolite detection rate was <80%, or were identified as outliers by multidimensional scaling analysis within a specific race/ethnic group. Metabolites were filtered if their detection rate across samples was <80% and, if applicable, had a coefficient of variation >20% for quality control (QC) samples. After quality filtering, missingness of each metabolite were imputed using the half minimum value, and the data were then standardized for analysis. Across all cohorts, we matched metabolites by their HMDB ID and/or PubChem ID, provided by the corresponding metabolomic laboratories. A total of 1,273 named metabolites were initially qualified for analysis in at least one cohort. To reduce single-study bias, we limited our analyses to 469 metabolites that were available in at least four independent cohorts, or available in at least three independent cohorts if the three cohorts covered both Metabolomic platforms. Finally, 407 metabolites from NHS, 363 from NHS2, 291 from HPFS, 364 from WHI, 327 from MESA, 274 from PREDIMED, 188 from FHS, 283 from SOL, 139 from ARIC and 231 from BPRHS were harmonized for our analysis (Extended Data Table 1). In CHS, 411 metabolites were included in genetic analyses (Supplementary Table 7). Details of the metabolomic profiling, QC and data processing are in the Supplementary Methods.
Metabolome-wide association analysis for incident T2D
Details of analytical approaches and models are provided in Supplementary Methods and Supplementary Table 1. Briefly, all association analyses were conducted separately for each cohort, stratified by major racial/ethnic groups when sample sizes permitted. Metabolites were inversely normal transformed by each substudy and racial/ethnic group (if applicable) in each cohort. To analyze the association between each metabolite and T2D risk, we applied Cox regression for studies of longitudinal cohort design (NHS excluding the T2D nested case–control substudy, NHS2, HPFS, SOL, ARIC, WHI, FHS, MESA and BPRHS); logistic regression for the NHS T2D nested case–control substudy; and Cox regression with Barlow weights80 and robust estimators for the PREDIMED T2D nested case–cohort study. The basic multivariate model (model 1) was adjusted for age, sex, smoking status, alcohol consumption and, if applicable, education, family income, fasting status, lipid-lowering medications, anti-hypertensive medications, family history of diabetes, self-reported physician-diagnosed hypertension, self-reported physician-diagnosed dyslipidemia and study-specific covariates. The main model was further adjusted for BMI and WHR (model 2). In sensitivity analyses, model 1 was further adjusted for PA and dietary quality index (model 3); high-density lipoprotein (HDL)-cholesterol, low-density lipoprotein (LDL)-cholesterol and triglycerides (model 4), or systolic and diastolic blood pressures (model 5). In another sensitivity analysis, model 2 was further adjusted for eGFR in NHS, NHS2, HPFS, SOL, ARIC, WHI and PREDIMED. For each metabolite, association results from all available cohorts and racial/ethnic groups were combined using a fixed-effect, inverse-variance-weighted (IVW) meta-analysis, and a meta-analyzed FDR < 0.05 was considered statistically significant. In secondary analyses, meta-analysis was conducted combining results from the same racial/ethnic groups, or cohorts using the same platforms.
To annotate the novelty of the identified associations, we reviewed previous prospective cohort studies linking circulating metabolites to T2D risk. We used a literature-review-based meta-analysis4 that included all studies published before 6 March 2021 as an anchor, and searched for additional studies published from 2021 to 202421,81,82,83,84,85,86,87,88,89,90,91,92,93,94. We considered an association as ‘previously reported,’ if the association was statistically significant in a published study after multiple testing correction based on the study’s prespecified analysis plan.
GWAS of metabolites
Detailed information on genotyping arrays, imputation methods, sample size and GWAS and meta-analysis methods, is provided in Supplementary Methods and Supplementary Table 7. Briefly, genotyping were conducted using several types of array by previous studies in NHS/HPFS95, SOL96, ARIC7, WHI97, CHS98 and FHS43. Imputation was conducted based on the HRC reference panel in NHS/HPFS and CHS; 1000 Genomes Project phase 3 worldwide reference panel in SOL, 1000 Genomes Project phase 3 v.5 in WHI and HapMap CEU population release v.22 in FHS with comprehensive pre- and postimputation QC. GWAS of metabolites were conducted previously in the NHS/HPFS (median n = 6,610, range 971–8,054) and WHI (n = 1,256) using the RVTESTS tool6,42,99, in SOL (n = 3,933) using a linear mixed-effect model in GMMAT7 and in ARIC (n = 1,772 and n = 1509 for African American and non-Hispanic white participants, respectively)7, CHS (n = 263) and FHS (n = 1,802)43, with detailed analysis procedures described in previous publications7,42,43.
GWAS summary statistics from each cohort were lifted over to Genome Build v.37 and filtered, retaining single nucleotide polymorphisms with a minor allele frequency ≥ 0.01 and imputation ratio ≥0.3. For each metabolite, an IVW fixed-effect meta-analysis, implemented in METAL100, was used to combine GWAS results from the cohorts in which the metabolite was available. Genomic control was implemented before and after meta-analysis100. The final GWAS were available for 458 out of 469 harmonized metabolites, with the total sample size ranging from 1,074 to 18,590 (median n = 8,611). We compared significant mQTLs identified at P < 5 × 10−8 and 1.09 × 10−10 (that is, 5 × 10−8 further correcting for 458 metabolites) levels. Manhattan plots were derived using R package CMplot and regional plots were draw with LocusZoom101. In a secondary analysis, we compared genetic effect heterogeneity between racial/ethnic groups at the identified mQTLs for T2D-associated metabolites (Supplementary Methods).
We annotate the novelty of our significant mQTLs for the 165 T2D-associated metabolites at P < 1.09 × 10−10, by comparing our results to eight previous studies (with N ≥ 4,000 and used LC–MS based metabolomic platforms)8,24,25,26,102,103,104,105. We considered a locus for a specific metabolite as ‘previously reported’ if the reported lead genetic variant was the same lead variant, or not the same lead variant but was significant in our study; or not in our study but within the clumping range of our identified locus. We considered a locus for a metabolite as potentially new if our locus was not previously reported for this metabolite, or this metabolite was not previously reported in these studies.
Lead variants for metabolites, pathway analysis and proportion of variance explained
We used the PLINK clumping function (P < 5 × 10−8 and r2 < 0.01 in a 1,000-kb window) to identify independent genetic variants associated with each metabolite. For metabolite with no variant at P < 5 × 10−8, a single lead variant with the smallest P was selected. Gene annotation for top variants was conducted using the SNPNexus web tool106. Canonical pathway enrichment analyses was conducted using the MetaCore software with the default background107; and we compared top enriched pathways for genes annotated to mQTLs of T2D-related metabolites versus those of non-associated metabolites. We calculated the R2 of each metabolite explained by independent lead genetic variants using the formula \({\sum }_{i=1}^{k}\beta \times \beta \times 2\times {\rm{MAF}}\times (1-{\rm{MAF}})\), in which k is the number of independent lead variants, and β is the association coefficient between the variant and the metabolite. We compared the R2 distribution for the T2D-associated versus non-associated metabolites using Wilcoxon test.
Genetic correlation r
g between metabolites and T2D-related traits
We acquired publicly available GWAS summary statistics from large consortium studies for T2D (180,834 cases and 1,159,055 controls)27, fasting insulin (N = 98,210)108, proinsulin (N = 45,861)109, HOMA-IR and HOMA-B (N = 51,750)110, BMI-adjusted insulin sensitivity index (ISI, N = 53,657) and insulin fold-change (IFC; N = 55,124)111, BMI and WHR (N = ∼700,000)112 and lipids (N = ∼1,500,000)113. We conducted GWAS for HBA1c (N = 390,982), subcutaneous fat volume (N = 37,912), visceral fat volume (N = 37,912), liver proton density fat fraction (PDFF; N = 29,512), pancreas PDFF (N = 28,624) and liver enzymes (N = ∼390,000) in the UK Biobank using BOLT-LMM (Supplementary Methods). We calculated rg between each metabolite and each clinical trait using linkage disequilibrium score regression, based on their GWAS summary data overlapping with the 1.2 M HapMap3 variants after excluding the major histocompatibility complex region in the European population114. For each trait, we compared the distribution of its rg with T2D-associated versus non-associated metabolites, using chi-squared test, and considered FDR < 0.05 (correcting for numbers of comparisons tested) as statistically significant.
Genetic colocalization
We obtained tissue-specific cis-eQTLs summary statistics from the GTEx project v.8115,116. The shared causal variants between each metabolite and tissue-specific transcriptome from 47 tissue types, were examined using colocalization analysis implemented in the coloc.abf() function in R package ‘coloc’ v.5117. For each metabolite, we input the GWAS summary statistics for all variants within ±500 kb of its independent lead variants (Supplementary Methods). A posterior probability of H4 (PPH4) > 0.8 was considered as strong evidence for genetic colocalization. Within each tissue type, we used univariant logistic regression to test whether the proportions of mQTL–eQTL colocalizations are higher for the T2D-associated versus non-associated metabolites, and a one-sided FDR < 0.05 (correcting for 47 tissue types) was considered as statistically significant. We applied a similar coloc approach to examine genetic colocalizations between circulating metabolites and T2D27. We then aligned mQTL–T2D colocalizations with tissue-specific eQTL–mQTL colocalizations by metabolites and shared causal variants, to interpret the potential functionality of metabolites in T2D pathogenesis.
MR analysis
To infer the potential causal relationships between 233 T2D-associated metabolites (with genetic data) and T2D risk, we applied four MR methods implemented in the MendelianRandomization R package118: we used mode-based estimate (MBE) as the main method as it is generally conservative and robust to outliers; we further applied weighted-median, IVW and MR-egger to indicate result consistency119. When testing the direction from metabolites to T2D, we used independent variants from clumping (P < 5 × 10−8 and r2 < 0.01 in a 1,000-kb window) excluding the HLA region as genetic instrumental variables. If fewer than three variants were identified, we reduced the clumping P threshold until at least three variants were identified. We considered a potential causal relationship when MBE–FDR < 0.05 and at least two other MR methods showed the same effect directions as those from MBE. Sensitivity analyses were conducted, either to remove variants mapped to the top 3 recurrent loci (GCKR, ZNF259, FADS cluster) from the instrumental variables, or to use only independent variants clumped at P < 1.09 × 10−10 as the instrumental variables of metabolites, using the IVW MR method (due to fewer variants retained). When testing the direction from T2D to metabolites, we used independent lead variants associated with T2D at P < 5 × 10−8 as the instrumental variables. For the 148 metabolites that are potential mediators between BMI and T2D risk, we applied MR analysis to test the direction from BMI to metabolites. Details are provided in Supplementary Methods.
MWASs for modifiable risk factors
We fitted linear models to regress inversely normal transformed metabolite levels on age, sex (only in SOL), current smoking status, BMI, PA, intakes of 15 main food groups and fasting status, simultaneously together with cohort-specific covariates. Analyses were conducted in NHS/HPFS, SOL and WHI, separately, further stratified by substudies or racial groups (Supplementary Methods). Association coefficients between metabolites and each particular risk factor were then combined across analytical sets using a fixed-effect IVW meta-analysis. The R2 of each metabolite explained by specific risk factors were first calculated in each analytical set using the formula \(\beta \times \beta \times {\mathrm{variance}}\left({\mathrm{risk}}\; {\mathrm{factor}}\right)\!{/\mathrm{variance}}\left({\mathrm{metabolite}}\right)\), with the β being the association coefficients between the metabolite and the risk factor; and then averaged across all analytical sets. We compared the distributions of R2 for T2D-associated versus non-associated metabolites using the Wilcoxon test.
Mediation analysis between risk factors, metabolites and T2D risk
Details for mediation analysis are described in Supplementary Methods. Briefly, our analysis focused on BMI, PA, coffee/tea consumption and red/processed meat intake. For each risk factor, metabolites (1) that were associated with both the risk factor and T2D risk and (2) whose association directions with the risk factor and T2D risk were consistent with the pre-assumed epidemiological relationships between the risk factor and T2D risk, were considered. We tested whether, and to what degree, each metabolite mediated the association between a risk factor and T2D risk using the CMAverse R package120, adjusting age, sex, smoking, BMI and PA (if not the tested risk factor), calorie intake and other cohort-specific covariates, separately in NHS/HPFS, SOL and WHI. We combined total, indirect and direct effects, respectively, from each analytical set using a fixed-effect meta-analysis. The mediated proportion was calculated by dividing indirect effect to total effect. Metabolites with an indirect effect FDR < 0.05 and a consistent effect direction between the indirect and total effects, was considered as a potential mediator between a risk factor and T2D risk.
A multimetabolite signature for incident T2D prediction
We used metabolites shared between the Broad Institute and the Metabolon platforms (excluding glucose) to develop the signature to increase its generalizability to future studies. To avoid overfitting in model development and testing, we employed a leave-one-cohort-out cross-validation approach, in which we set aside one cohort as the testing set each time, and trained a prediction model for the set-aside cohort using data from all other cohorts (Extended Data Fig. 8). Given the heterogeneity of our cohorts, we did not pool individual-level data for model training. Instead, we applied a two-step approach to train the prediction model in a representable cohort (that is, WHI, which assessed the most shared metabolites for all its participants) but also leveraged association data from several other cohorts. In each iteration (that is, for each held-out testing cohort), we first conducted a metabolome-wide meta-analysis for T2D risk using all cohorts except WHI and the held-out cohort. Then, metabolites associated with T2D risk at FDR < 0.05 in the first step and shared between the two metabolomic platforms, were used as input in a Cox regression with elastic net regularization, implemented using the glmnet R package121, to construct a metabolomic signature model for T2D prediction in WHI. The derived model was further applied to the held-out cohort to calculate a metabolomic signature score. Within WHI, a leave-one-out cross-validation approach was used to acquire the unbiased metabolomic signature score. For details, please see Supplementary Methods.
The metabolomic signature scores, calculated in each held-out cohort, were then standardized. To evaluate whether the signature improved the T2D risk prediction, we fitted three sets of logistic (in SOL, and T2D nested case–control substudy in NHS) or Cox models (all other datasets): one model including only the metabolomic signature; a conventional risk factor model including age, sex, smoking, lipid-lowering medication use, anti-hypertensive medications, family history of diabetes, hypertension, dyslipidemia and BMI; and a third model including all conventional risk factors and the metabolomic signature. We compared the AUC between the conventional model versus the conventional plus metabolomic signature model. In a secondary analysis, we further included blood glucose (from metabolomic assays) in the conventional model to evaluate the added value of the metabolomic signatures beyond blood glucose.
In each cohort, we calculated the crude incident rate of T2D across deciles of the signature score. We fitted logistic or Cox models to analyze the relative risk of T2D, comparing higher versus lowest deciles of the metabolomic signature, adjusting for the same covariates in the main analysis model 2. In NHS/HPFS, SOL and WHI, we examined associations between the metabolomic signature with baseline risk factors, by regressing the signature score on age, sex (if appropriate), current smoking status, BMI, PA, intakes of 15 main food groups and fasting status simultaneously, together with cohort-specific covariates, using linear regression. All analysis was conducted separately in each cohort, and results were combined using a meta-analysis. FDR < 0.05 was considered as statistically significant.
We conducted two sensitivity analyses during model development. One was to use SOL (measured the most metabolites using the Metabolon platform) as the representative training cohort instead of WHI, which showed a similar, albeit slightly weaker, model performance in held-out cohorts (Extended Data Fig. 8). The other was to compare between elastic net versus lasso regularizations121, which reaffirmed that elastic net regression had compatible but a slightly better performance versus lasso regression (Supplementary Fig. 13). Separately from the leave-one-cohort-out cross-validation, we presented a final metabolomic signature model for future studies, developed using data from all study cohorts. For this model, we first conducted a metabolome-wide meta-analysis for T2D risk in all cohorts except WHI, and then used significant metabolites (FDR < 0.05) as input in a Cox regression with elastic net regularization for T2D prediction in WHI. The selected metabolites and their coefficients of this final model are highly consistent with those of models applied to each held-out cohort (Supplementary Table 18a).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.