Multivariate genome-wide analysis of aging-related traits identifies novel loci and new drug targets for healthy aging

Abstract

The concept of aging is complex, including many related phenotypes such as healthspan, lifespan, extreme longevity, frailty and epigenetic aging, suggesting shared biological underpinnings; however, aging-related endpoints have been primarily assessed individually. Using data from these traits and multivariate genome-wide association study methods, we modeled their underlying genetic factor ('mvAge'). mvAge (effective n = ~1.9 million participants of European ancestry) identified 52 independent variants in 38 genomic loci. Twenty variants were novel (not reported in input genome-wide association studies). Transcriptomic imputation identified age-relevant genes, including VEGFA and PHB1. Drug-target Mendelian randomization with metformin target genes showed a beneficial impact on mvAge (P value = 8.41 × 10^-5). Similarly, genetically proxied thiazolidinediones (P value = 3.50 × 10^-10), proprotein convertase subtilisin/kexin 9 inhibition (P value = 1.62 × 10^-6), angiopoietin-like protein 4, beta blockers and calcium channel blockers also had beneficial Mendelian randomization estimates. Extending the drug-target Mendelian randomization framework to 3,947 protein-coding genes prioritized 122 targets. Together, these findings will inform future studies aimed at improving healthy aging.

Fig. 1. Study overview.

An overview of this study’s data sources, analytical flow and methodology. Created with BioRender.com. The univariate input GWASs of frailty and EAA were reverse coded to align their effects to have positive relationships with healthspan, lifespan and extreme longevity. GWAS, genome-wide association study; EAA, epigenetic age acceleration; CELLECT, CELL-type Expression-specific integration for Complex Traits; IVW, inverse variance weighted; MR LASSO, MR Least Absolute Shrinkage and Selection Operator; HbA1c, glycated hemoglobin; LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol.

Fig. 2. Multivariate aging GWAS modeled with genomic SEM.

a, Genetic correlations for SEM with genomic SEM, displaying pairwise LD score genetic correlation estimates for the five univariate phenotypes. b, Path diagram of the common factor model estimated with genomic SEM, with standardized factor loadings (standard error in parentheses). c, Manhattan plot showing SNP associations (−log₁₀(P value)) with mvAge, ordered by chromosome. The red dashed line indicates the threshold for conventional genome-wide significance (P value = 5 × 10⁻⁸). P values are derived from two-sided Wald tests for each SNP on mvAge. * indicates that summary statistics for frailty and PhenoAge (the epigenetic clock variable) were reversed to align with the other longevity-related endpoints. µ reflects the residual variance in the genetic indicators for the input univariate age-related GWASs not explained by the mvAge common factor. Source data

Fig. 3. Drug-target MR results assessing proxying metformin and other antidiabetics targets with mvAge.

Data presented are MR effect estimates (betas) for the IVW MR method (the primary MR method) and the corresponding 95% confidence intervals (CIs) aligned to proxy the pharmacological effect of metformin and antidiabetic genes (HbA1c GWAS n = 344,182) (1 s.d. lowering of HbA1c levels) on mvAge (n = 1,958,774). The vertical line in the center of the forest plots is 0, corresponding to no change in the IVW estimate of the drug targets on mvAge. Full results are presented in Supplementary Tables 15 and 16. Metformin results plotted show the MR estimates for the primary metformin instrument (top row), which comprised variants within genes for five metformin targets (AMPK, GSD15, MCI, MG3 and GLP1), for the estimates of the alternative metformin instruments used as sensitivity analyses, and the metformin targets separated into individual instruments. See Methods and Supplementary Methods sections for additional details. For the analyses of the antidiabetic classes not including the metformin targets, the * indicates that the thiazolidinedione MR estimate surpasses the Bonferroni-adjusted P-value threshold = 0.002, corrected for the 25 antidiabetic, lipid-modulating and antihypertensive drug targets compared. P values are derived from two-sided Wald tests. Gene names for the nearest mapped genes are italicized. Source data

Fig. 4. Drug-target MR results assessing the impact of lipid-modulating and antihypertensive gene targets on mvAge.

Data presented are MR effect estimates (betas) for the IVW MR method (the primary MR method) and the corresponding 95% CIs aligned to proxy the pharmacological effect of modulated lipid levels (1 s.d. lower LDL-C (n = 440,546), 1 s.d. lower TG (n = 441,016), and 1 s.d. higher HDL-C (n = 403,943)) and SBP (n = 436,419) (antihypertensive gene targets (1 s.d. lower SBP) on mvAge (n = 1,958,774)). The vertical line in the center of the forest plots is 0, corresponding to no change in the IVW estimate of the drug targets on mvAge. Full results are presented in Supplementary Table 17. * indicates a P value surpassing the Bonferroni-adjusted P-value threshold = 0.002, corrected for the 25 antidiabetic, lipid-modulating and antihypertensive drug targets compared. P values are derived from two-sided Wald tests. Gene names for the nearest mapped genes are italicized. Lp(a), lipoprotein a. Source data

Fig. 5. Cis-instrument MR results assessing the impact of protein-coding genes on mvAge through their associations with HbA1c and LDL-C.

a, Volcano plot of the Z scores (versus the negative log₁₀(P value)) of the MR estimates (beta/se) for the inverse variance weighted MR method aligned to proxy the pharmacological effect of lowered HbA1c levels. b, Volcano plot of the Z scores (versus the negative log₁₀(P value)) of the MR estimates (beta/se) for the inverse variance weighted MR method aligned to proxy the pharmacological effect of lowered LDL-C levels. Dotted lines indicate the Bonferroni-corrected P-value threshold (1.92 × 10⁻³). Labeled genes are those with beneficial estimates on mvAge that surpass the Bonferroni-corrected P-value threshold and align with lowered HbA1c and lower LDL-C. c, The STITCH protein–protein and protein–chemical interactions for the 30 protein-coding genes in HbA1c. Stronger associations are annotated with thicker lines. Protein–protein interactions are represented by gray lines, protein–chemical interactions are represented by green lines, and chemical–chemical interactions are represented by red lines. d, Flowchart outlining the cis-instrument analysis pipeline (see Supplementary Methods for more details). P values are derived from two-sided Wald tests. Source data

Extended Data Fig. 1. Gene-level Manhattan plot of transcriptome-wide association study.

TWAS Z-scores (effect estimates of the imputed gene expression on mvAge) are plotted (mvAge N = 1,958,774). The expression quantitative trait loci (eQTL) data are derived from the GTEx Version 8 sparse canonical correlation analysis performed by Feng et al. 2021 (see References). The blue lines represent Z = ± 4.69, corresponding to the Bonferroni adjusted P-value threshold = 1.38 × 10⁻⁶ (0.05/36,149 sCCA features available for analysis). Red points and labels indicate genes (Ensembl gene IDs) surpassing the threshold. TWAS, transcriptome-wide association study; sCCA, sparse canonical correlation analysis. Source data

Extended Data Fig. 2. Cis-instrument Mendelian randomization assessing the impact of protein coding genes on mvAge through their associations with HDL-C, LDL-C, triglycerides, and SBP.

Presented is a volcano plot of the Z scores (versus the negative log₁₀(P-value)) of the MR estimates (beta/se) for the inverse variance weighted (IVW) MR method aligned to proxy the pharmacological effect of increased HDL-C (N = 403,943), lowered LDL-C (N = 440,546), lowered TG (N = 441,016), or lowered SBP (N = 436,419) on mvAge (N = 1,958,774). The dotted line indicates the Bonferroni corrected P-value threshold (1.92 × 10⁻³). Labeled genes are those with beneficial estimates on mvAge that surpass Bonferroni corrected P-value threshold and align with increased HDL-C, lowered LDL-C, lowered TG, or lowered SBP. HDL-C, high-density lipoprotein cholesterol; LDL-C, low-density lipoprotein cholesterol; TG, triglycerides; SBP, systolic blood pressure; MR, Mendelian randomization. Source data

Extended Data Fig. 3. Protein-protein and protein-chemical interaction for protein coding genes located near the lead SNPs of the HDL-C, LDL-C, triglycerides, and SBP GWASs associated with mvAge in the cis-instrument Mendelian randomization.

Analyses were performed using STITCH (http://stitch.embl.de/). Results plotted are STITCH PPI interaction scores. The combined scores are computed by calculating the probabilities from the STITCH database sources of evidence and correcting by the probabilities of randomly observing interactions between the proteins and chemicals (see ref. in the Reference list for additional information). Stronger associations are annotated with thicker lines. Protein-protein interactions are represented by grey lines, protein-chemical interactions are represented by green lines, and chemical-chemical interactions are represented by red lines. Per STITCH guidelines, ‘high confidence’ interactions are considered those with combined scores ≥ 0.90. HDL-C, high-density lipoprotein cholesterol; LDL-C, low-density lipoprotein cholesterol; TG, triglycerides; SBP, systolic blood pressure; MR, Mendelian randomization. Source data

Extended Data Fig. 4. Cis-instrument Mendelian randomization analysis of 68 circulating proteins.

Data presented are MR effect estimates (betas) for the inverse variance weighted (IVW) MR method and corresponding 95% confidence intervals (CIs). The impact of 68 circulating proteins on mvAge (N = 1,958,774) was analyzed using protein quantitative trait loci (pQTL) data derived from 30,391 participants of European ancestry in the SCALLOP Consortia dataset (http://www.olink-improve.com/). We used pQTLs associated with the respective plasma protein at P-value < 5 × 10⁻⁸ within or near the cis-acting locus of the target gene boundary, that is, 100 kilobases on either side of the respective gene boundary. Extracted SNPs were clumped at the LD R² ≤ 0.2 threshold (250 kb) using the 1000 Genomes Phase 3 European reference population and MR IVW (random-effects analysis when there were more than three variants) and MR Egger implemented accounting for correlation between the instrument variants. MR, Mendelian randomization. Source data

Extended Data Fig. 5. Mendelian randomization model overview (directed acyclic graph).

B₂ is the genetic association of interest, estimated by B₂ = B₁/ B_3. B₁ and B₃ represent the estimated MR association of the genetic variants on the exposure and the outcomes. MR assumes that the genetic variants comprising the instrument for the exposure only impact the outcome of interest via the exposure and not directly, or via confounders (dotted lines). MR, Mendelian randomization.

Extended Data Fig. 6. Drug-target Mendelian randomization analysis overview of anti-diabetics, lipid-modulating targets, and antihypertensives.

Diagram depicts flow diagram and details of the drug-target MR analyses of the cardiometabolic targets on mvAge (N = 1,958,774). Prior to step 1, colocalization analysis was employed to prioritize protein-coding targets for the screen. See Methods and Supplementary Methods for details on selection and identification of the individual targets in the three broad drug classes (anti-diabetics, lipid-modulating targets, and antihypertensive). In step 2, cis-instrumentation was performed using genome-wide association study (GWAS) of biomarkers that are the primary indications of pharmacological modulation of these targets. For antidiabetics, GWAS data of HbA1c (N = 344,182) was used; for lipid-modulating targets, several lipid subfractions including LDL-C (N = 440,546), triglycerides (N = 441,016), and HDL-C (N = 403,943) were used; and for antihypertensives, GWAS data of SBP (N = 436,419). Independent variants (LD R² < 0.2) at P-values < 5 × 10⁻⁸) were extracted, and cis-instruments constructed for each target, which exposure variants were then extracted from the mvAge GWAS (outcome), harmonized, and then analyzed using multiple MR methods (steps 3 and 4). See Methods and Supplementary Methods for further details. MR, Mendelian randomization; LD, linkage disequilibrium. LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol; TG, triglycerides; SBP, systolic blood pressure.

Extended Data Fig. 7. Analysis overview of cis-instrument Mendelian randomization screen of protein coding genes associated with either HbA1c, LDL-C, HDLC, triglycerides, or SBP.

Results show number of protein-coding genes including in the stages of the colocalization and cis-instrument MR screens. Exposure genome-wide association study data for each of the biomarkers came from HbA1c (N = 344,182); several lipid subfractions including LDL-C (N = 440,546), triglycerides (N = 441,016), and HDL-C (N = 403,943); and SBP (N = 436,419). Outcome data was mvAge (N = 1,958,774). HbA1c, glycated hemoglobin; HDL-C, high-density lipoprotein cholesterol; LDL-C, low-density lipoprotein cholesterol; TG, triglycerides; SBP, systolic blood pressure; MR, Mendelian randomization.