Mendelian randomization of genetically independent aging phenotypes identifies LPA and VCAM1 as biological targets for human aging

Length and quality of life are important to us all, yet identification of promising drug targets for human aging using genetics has had limited success. In the present study, we combine six European-ancestry genome-wide association studies of human aging traits—healthspan, father and mother lifespan, exceptional longevity, frailty index and self-rated health—in a principal component framework that maximizes their shared genetic architecture. The first principal component (aging-GIP1) captures both length of life and indices of mental and physical wellbeing. We identify 27 genomic regions associated with aging-GIP1, and provide additional, independent evidence for an effect on human aging for loci near HTT and MAML3 using a study of Finnish and Japanese survival. Using proteome-wide, two-sample, Mendelian randomization and colocalization, we provide robust evidence for a detrimental effect of blood levels of apolipoprotein(a) and vascular cell adhesion molecule 1 on aging-GIP1. Together, our results demonstrate that combining multiple aging traits using genetic principal components enhances the power to detect biological targets for human aging. Many aging-related phenotypes share a common genetic component, but to disentangle disease-specific variants from aging-specific ones has been challenging. Here Timmers et al. combined several genetics studies of aging-related traits to identify common underlying genetic factors that contribute to aging.

A ging affects us all, from the personal, progressive loss of health to the collective burden of chronic age-related disease and frailty on society. In humans, the body undergoes a systemic functional decline after reaching adulthood, which manifests itself as age-related disease, infirmity and eventually death 1 . The factors determining the rate of aging and death are complex and interlinked, and include genetics, lifestyle, environmental exposures and chance.
Quantifying the aging process is not straightforward. A variety of aging-related phenotypes has been studied as proxies, from chronological measurements such as the length of time from birth until the occurrence of a major disease (healthspan) 2 or death (lifespan) 3,4 , to cellular deterioration measurements such as telomere attrition 5 and loss of the Y chromosome 6,7 , to holistic measurements such as the frailty index 8 , encompassing multiple functional impairment indicators 9 . Although the genetic component of these aging-related traits tends to be estimated at <15% 2,8,10 , recent progress has been made on characterizing this component using large genome-wide association studies (GWASs), and combining similar GWASs to increase statistical power 11,12 .
The benefit of combining GWASs of several aging phenotypes, especially in different populations, is the ability to detect biological mechanisms that influence multiple core components of aging, while downweighing population-and trait-specific features. For example, a recent multivariate analysis of healthspan, parental lifespan and longevity GWASs found that genetic loci that were not shared between traits often associated with population-specific, behavioral risk factors such as smoking and skin cancer 11 . On the other hand, that analysis showed that genetic loci shared between traits were associated with biological pathways such as cellular homeostasis and heme metabolism 11 .
However, to date, large aging-related trait GWASs have been combined only using multivariate analysis of variance, which detects Mendelian randomization of genetically independent aging phenotypes identifies LPA and VCAM1 as biological targets for human aging Length and quality of life are important to us all, yet identification of promising drug targets for human aging using genetics has had limited success. In the present study, we combine six European-ancestry genome-wide association studies of human aging traits-healthspan, father and mother lifespan, exceptional longevity, frailty index and self-rated health-in a principal component framework that maximizes their shared genetic architecture. The first principal component (aging-GIP1) captures both length of life and indices of mental and physical wellbeing. We identify 27 genomic regions associated with aging-GIP1, and provide additional, independent evidence for an effect on human aging for loci near HTT and MAML3 using a study of Finnish and Japanese survival. Using proteome-wide, two-sample, Mendelian randomization and colocalization, we provide robust evidence for a detrimental effect of blood levels of apolipoprotein(a) and vascular cell adhesion molecule 1 on aging-GIP1. Together, our results demonstrate that combining multiple aging traits using genetic principal components enhances the power to detect biological targets for human aging.
tion matrix (Fig. 2), yielding association summary statistics for six aging-GIPs (available at https://doi.org/10.7488/ds/2972). As expected, high-definition likelihood (HDL) 21 estimated aging-GIP1 to be the most heritable of the aging-GIPs (h 2 SNP = 0.20; standard error (s.e.) = 0.005), capturing >70% of the genetics of healthspan, parental lifespan and longevity, and 90% of the genetics of selfrated health and the inverse of frailty (henceforth referred to as 'resilience') ( Fig. 3). A leave-one-out analysis confirms that the GIP analysis is highly robust to the selection of GWAS, with the genetic architecture of aging-GIP1 remaining largely the same, after excluding any one of the six chronological and holistic aging-related trait GWASs (estimates range from r g = 0.950 (s.e. = 0.024) when excluding resilience, to r g = 0.996 (0.023) when excluding healthspan; Supplementary Table 2).
Apart from the genetic correlations with its component traits, aging-GIP1 shows significant correlations, with 459 of 729 traits tested (P adj < 0.05, Bonferroni-adjusted for 6 GIPs and 729 traits). These suggest that, in addition to length of life, aging-GIP1 also captures mental and physical wellbeing. For example, among the largest negative aging-GIP1 correlations are mental illness, taking medications and diseases of old age (such as cardiometabolic disorders, cancers and osteoarthritis; r g ≤ -0.5; P adj < 1 × 10 −9 ). Inversely, fitness and education are among the largest positive correlations (r g ≥ 0.5; P adj < 1 × 10 −10 ). Aging-GIP1 also shows moderate negative correlations with infectious diseases, including N39 (International Classification of Diseases, 10th revision 22 ) urinary tract infections (r g = -0.66; 95% CI = -0.42 to -0.90; P adj = 3 × 10 −4 ), coughing on most days (r g = -0.39; 95% CI = -0.30 L o n g e v i t y F a t h e r l i f e s p a n M o t h e r l i f e s p a n H e a l t h s p a n F r a i l t y H e a l t h ( s e l f -r a t e d ) to -0.49; P adj = 1 × 10 −11 ) and severe COVID-19 hospitalization (that is, resulting in respiratory support or death; r g = -0.33; 95% CI = -0.20 to -0.46; P adj = 0.004).
However, aging-GIP1 also retains some strong correlations with socioeconomic factors, such as smoking behaviors (for example, current tobacco smoking: r g = -0.50; 95% CI = -0.45 to -0.55;  ) and SNP heritability (h 2 SNP ) of each aging-GIP GWAS, with s.e. in parentheses. b, Genetic correlations between aging-GIP GWASs and the six aging-related trait GWASs used to construct them. Blank values failed to pass nominal significance (P ≥ 0.05). c, GWASs of 402 phenotypes measured in European-ancestry individuals from the GWAS-MAP platform with strong and significant associations with GIP1 (P < 1.5 × 10 −5 and |r g | > 0.25). These 402 traits are summarized here in 25 groups, based on hierarchical clustering of the magnitude of pairwise genetic correlations. Each group is manually annotated with a label describing the traits within the cluster, and displays the values of the most informative trait (that is, the trait with the highest total z-score across aging-GIPs). Values failing to pass nominal significance (P ≥ 0.05) are grayed out. The dendrogram displayed here shows the hierarchical relationship between the most informative trait in each cluster. HDL, high-density lipoprotein; VLDL, very-low-density lipoprotein. See Supplementary Data 3 for the full list of correlations. Genetic correlation P values are derived from two-sided Wald's test with one degree of freedom. P adj = 4 × 10 −92 ) and having a job involving manual or physical work (r g = -0.49; 95% CI = -0.44 to -0.54; P adj = 2 × 10 −89 ; Fig. 3 and Supplementary Data 3). Although such socioeconomic risk factors are partly heritable 23 , they could also reflect social stratification or geographic confounding shared between the original six studies used to construct aging-GIP1 (refs. 24,25 ). Therefore, as a sensitivity analysis, we used the GIP methodology to remove all genetic correlations with UK Biobank GWASs of household income and socioeconomic deprivation from the aging-GIPs-recognizing that this may remove some of the true signal as well (Extended Data Fig. 2). The SNP heritability of the adjusted aging-GIP1 phenotype is substantially attenuated (h 2 SNP = 0.09; s.e. = 0.003) and, consequently, genetic correlations with all phenotypes are reduced. However, most traits that were significantly associated at 5% FDR remain associated with the residual phenotype (Extended Data Fig. 3). As expected, the largest reductions in genetic correlations are with socioeconomic, educational and smoking/drinking traits, which are reduced by ≥0.30 (Supplementary Data 3).
The remaining aging-GIPs have much lower heritability (h 2 SNP < 0.065) and are of less interest to the present study because they capture genetic variance shared between fewer traits and do so with higher degrees of uncertainty. In short, aging-GIP2 correlates mainly with nonlethal causes of poor self-rated health, such as chronic pain (neck, back, stomach), negative emotions and hearing problems (59 traits; P adj < 0.05); aging-GIP3 correlates mainly with measures of socioeconomic deprivation and body composition (73 traits); and aging-GIP4 captures healthspan-specific correlations with cancer and diabetes (7 traits). Although aging-GIP5 and aging-GIP6 are underpowered due to their low heritabilities, they appear to distinguish parental lifespan from longevity (possibly through educational attainment), and between father lifespan and mother lifespan (possibly through risk taking and cardiovascular factors), respectively (Fig. 3).
Characterizing the genomics of aging-GIP1. Across the genome, 27 loci in the aging-GIP1 GWAS pass the genome-wide significance threshold, Bonferroni's adjusted for 6 GIPs (P < 8.3 × 10 −9 ). The strongest lead SNPs in these loci are rs429358 (P = 2 × 10 −40 ) and rs660895 (P = 2 × 10 −22 ), located nearest to the APOE and HLA-DRB1/DQA1 genes, respectively ( Fig. 4 and Table 1). Genetic fine-mapping with SuSiE 26 highlights three independent 95% credible sets of causal variants for the APOE locus: the APOE ε4 and ε2 alleles at single variant resolution, and a set of three variants intronic or near to APOC1. Three other loci have a credible set containing only one variant. These variants are rs9271300 (intergenic between HLA-DRB1 and HLA-DQA1), rs34811474 (a nonsynonymous ANAPC4 exon variant) and rs2165702 (an intergenic variant near PHB; Supplementary Table 3 and Supplementary Data 4).
Most lead SNPs are strongly associated with self-rated health and resilience, in line with the large loadings of these traits in the construction of aging-GIP1 ( Fig. 4 and Supplementary Table 4). Seventeen aging-GIP1 loci overlap with genome-wide significant (P < 5 × 10 −8 ) loci from the original aging-related trait GWAS, showing evidence of a shared causal variant (coloc posterior probability (PP) ≥ 80%; Supplementary Table 5). For the other loci, the GIP framework either increases power (n = 4) or indicates that there may be a different causal variant (coloc PP < 80%; n = 6; Supplementary Data 5).
Loci near APOE, HLA-DRB1/DQA1, LPA and CDKN2B/-AS1 have previously been validated using the same trait in an external cohort [2][3][4] . For the remaining 23 loci, we measured lead SNP effects on participant survival in FinnGen (release 5; n = 203,244; 6.94% deceased) and BioBank Japan (n = 135,983; 24.1% deceased), to provide additional evidence of their association with human aging traits in independent samples. Combining both cohorts to achieve adequate power, we find that aging-GIP1-increasing alleles of lead SNPs near HTT and MAML3 have a protective effect on survival in these cohorts (one-sided P < 0.0022, that is, taking into account the 23 loci tested), increasing average lifespan by around 2.53 (95% CI = 0.91-4.15; one-sided P = 0.001) and 2.51 (95% CI = 0.79-4.23; one-sided P = 0.002) months per allele, respectively. Although we are underpowered to confirm the remaining 21 loci individually, we find that, collectively, their aging-GIP1-increasing alleles are also associated with increased Finnish and Japanese survival (one-sided P = 0.044; Supplementary Table 6).
Again, we find aging-GIP1 genetics are stable when performing a leave-one-out analysis: lead SNP effects of most aging-GIP1 loci do not change significantly when excluding one of the six aging traits. The exceptions are for previously replicated loci and SLC22A1/A2, where exclusion of one of the traits can reduce the aging-GIP1 effect size, although loci near HLA-DRB1/DQA1, LPA and CDKN2B/-AS1 remain genome-wide significant in all leaveone-out scenarios (Supplementary Data 6). Similarly, when completely removing genetic correlations with household income and socioeconomic deprivation, we find that aging-GIP1 effect sizes are attenuated but, in all cases, remain associated with the residual trait (P < 0.0019, that is, taking into account the 27 loci tested; Extended Data Fig. 4 and Supplementary Data 7).
We further looked up all lead aging-GIP1 SNPs and close proxies (r 2 EUR ≥ 0.8) in PhenoScanner 27 and the GWAS catalog 28 , excluding associations discovered solely in UK Biobank, which showed that 24 of 27 aging-GIP1 loci had previously been associated with one or more traits at genome-wide significance. Most of these loci were associated with cardiometabolic, immune-related or neuropsychiatric disorders, although several were also associated with measures of educational attainment and household income. Of specific interest are loci near APOE, HLA-DRB1/DQA1, CDKN2B/-AS1 and ZNF652/PHB, which show aging-GIP1-increasing alleles are associated with a reduction in multiple diseases, but do not appear to associate with socioeconomic factors, suggesting that these loci largely capture intrinsic sources of aging (Supplementary Data 8).
Aggregating SNP association statistics across the genome into gene scores using the Pathway Scoring Algorithm (PASCAL) 29 , we find that high-scoring genes for aging-GIP1 (Supplementary Data 9) appear to be overrepresented in the heme metabolism Hallmark gene set, as well as 300 gene ontology (GO) pathways (FDR = 5%), regardless of adjustment for socioeconomic correlations. These GO pathways cluster into 20 groups, related to: neuronal development, organization and function; transcriptional regulation; chemical homeostasis; cellular growth, differentiation and apoptosis; proteolysis; protein phosphorylation; intracellular signaling and transport; immune system development; the muscle system; and lipoprotein metabolism (Supplementary Data 10). Similarly, aging-GIP1 heritability appears to be enriched in genomic regions containing histone marks associated with the central nervous system (Supplementary Table 7).

Causal inference of blood protein levels on aging-GIP1.
Mendelian randomization (MR) uses genetic variation as instrumental variables to estimate the casual effect of exposures of interest on an outcome 30 . We used a set of well-validated blood protein quantitative trait loci (pQTLs) for 857 proteins 30 as genetic instruments in a two-sample MR 31 and colocalization 32 framework to infer putative causal links between protein levels and aging-GIP1 (Supplementary Data 11). We find robust evidence for a detrimental effect on aging-GIP1 (FDR < 5%) for the levels of four proteins in blood (Table 2), with pQTL instruments passing sensitivity and causal directionality tests, and additionally colocalizing with the aging-GIP1 signal (Methods). Three of these proteins-apolipoprotein(a) (LPA), olfactomedin-1 (OLFM1) and low-density lipoprotein (LDL) receptor-related protein 12 (LRP12)-were instrumented by a cis-pQTL and encoded by genes that appeared significantly enriched in the gene score analysis. The remaining protein, vascular cell adhesion molecule 1 (VCAM1), was instrumented by a trans-pQTL shared with β 2 -microglobulin (β 2 M); however, only VCAM1 protein levels colocalized with the signal at this locus (Supplementary Data 11). When performing the same MR analysis on aging-GIP1 adjusted for household income and socioeconomic deprivation, protein effects remained significant (Extended Data Fig. 5).
Among the discovered proteins, LPA shows the most significant effect on aging-GIP1 (P MR = 2 × 10 −8 ), with an increase of 1 s.d. in genetically predicted blood protein levels, causing a decrease of 0.035 (95% CI = 0.025-0.045) s.d. in aging-GIP1. This significance appears to be driven by a consistent detrimental effect across all six aging-GIP1 component traits (β MR range 0.013-0.035; all nominal P < 0.05; Fig. 5). For a sense of scale, when performing the same MR analysis on the unstandardized parental lifespan GWAS, this equates to a loss of approximately 7 months of life (95% CI = 5-9 months) per s.d. increase in LPA blood levels (Supplementary Table 8). The effects of the remaining proteins on aging-GIP1 are larger in magnitude but appear unequally distributed across aging-GIP1 component traits. We estimate that a decrease of 1 s.d. in genetically predicted VCAM1 blood protein levels is associated with an increase of approximately 18 months of life (13-23 months); however, comparing VCAM1 effects on standardized aging-GIP1 components, we find its effect to be largest on mid-to-late-life aging traits (lifespan, longevity, resilience) compared with early life aging traits (healthspan). Similarly, the effects of OLFM1 and LRP12 proteins appear mostly mediated through late-life aging traits (resilience and/or longevity; Fig. 5).

Discussion
We combined European-ancestry GWASs of healthspan, father lifespan, mother lifespan, longevity, frailty and self-rated health using a framework that maximized power to detect associations with their shared genetic component while downweighing trait-specific genetic associations. The resulting aging-GIP1 trait captured the genetics underlying physical and mental wellbeing, and showed strong inverse genetic correlations with cardiovascular, inflammatory and neuropsychiatric disease traits. We highlight 27 loci with genome-wide significant effects on aging-GIP1, including two new loci near HTT and MAML3 which showed directionally consistent evidence of an effect on survival in two additional, independent samples. Across the genome, we found genes enriched for association with aging-GIP1 to be overrepresented in heme metabolism and pathways related to (among others) neurogenesis, homeostasis, Resilience  proteolysis, immunity and the muscle system in human aging. Last, we performed MR of predicted blood protein levels on aging-GIP1, which revealed that the levels of LPA, VCAM1, OLFM1 and LRP12 proteins may be detrimental to multiple indices of healthy aging. LPA is a glycoprotein making up the main component of large lipoprotein(a) particles. It is a well-known risk factor for atherosclerotic disease 33,34 and a target of ongoing clinical trials with regard to cardiovascular outcomes 35 . A recent MR study used 27 genetic instruments for LPA and found a link between genetically elevated LPA levels and reduced healthspan and parental lifespan 36 , in line with our results. Our analysis suggested that the detrimental effect of LPA may apply to aging more generally, and stringent colocalization and reverse MR tests provided additional evidence for causality. It is interesting that human endothelial cell culture experiments Loci were defined as 500-kb regions centered on a lead genome-wide significant SNP (P < 8.3 × 10 −9 ) in linkage equilibrium (r 2 EUR < 0.1) with other lead locus SNPs. Nearest gene(s), closest genes upstream/ downstream to the lead SNP (within 250 kb) or, if none, the closest cytogenetic band; rsID, the lead SNP within the locus; Position, base-pair position (GRCh37); A1, effect allele, associated with higher aging-GIP1; Freq1, allele frequency of the effect allele in UK Biobank; β1, effect estimate (and s.e.) of the A1 allele on aging-GIP1 in s.d. units; P, two-sided nominal P value from Wald's test; Het, evidence of heterogeneity: asterisks indicate that aging-GIP1 effect size changes significantly when leaving out one of the core aging traits from the GIP1 calculation (*all but one leave-one-out effects are the same; **all but two effects are the same; ***≤3 effects are the same). show that the addition of LPA increases cell-surface expression of VCAM1 (ref. 37 ) and can increase endothelial cell contraction and permeability 38 . VCAM1 is a cell adhesion glycoprotein localized predominantly on endothelial cell surfaces and its expression is upregulated in response to inflammatory signals, which mediate adhesion and transduction of leukocytes across endothelial walls 39 . The link that we established between VCAM1 and human aging relied on a trans-pQTL instrument shared with β 2 M and is therefore more susceptible to horizontal pleiotropy, that is, the genetic variant may influence VCAM1 levels indirectly, and its effect on human aging traits could be caused by factors independent of VCAM1 levels. However, only VCAM1 colocalized with the pQTL signal, and experimental evidence from mouse studies suggests that the effect of VCAM1 on aging is likely to be causal. Specifically, VCAM1 levels in blood are known to increase with age in both humans and mice 40 , and treatment with anti-VCAM1 antibodies or an inducible deletion of Vcam1 improves cognitive performance of aged mice 40 . Of note, similar results have been found for β 2 M abundance and mouse knockouts 41 and, as such, identification of robust genetic instruments for β 2 M levels is also warranted.
We were unable to robustly assess colocalization and reverse causality for the LRP12 and OLFM1 signals because genome-wide association summary statistics were not available, so the effects of these proteins on human aging should be interpreted with additional caution. LRP12 belongs to the LDL-receptor superfamily and may play a role in brain development 42 and both tumor proliferation and suppression, depending on the tissue 43,44 . Similarly, OLFM1 is a glycoprotein involved in neuronal development and maintenance 45 , which has also been shown to suppress colorectal tumor metastasis 46 . In mouse models, Olfm1 knockouts showed reduced cerebral infarction and fertility 47 . However, it is unclear whether reduction of LRP12 or OLFM1 blood levels will have a beneficial effect in humans.
Importantly, our MR analysis was restricted to blood pQTLs, precluding detection of causal effects of protein levels in other tissues. Although it is likely that there are proteins with tissue-specific effects on aging, particularly in the brain, the samples needed for detection of such pQTLs are less readily available and therefore more difficult to study at scale (although progress is being made 48 ). Still, blood may be a particularly suitable tissue to identify agingrelated proteins. Connecting the circulatory systems of two mice of different ages has been shown to accelerate signs of aging in the brain, muscle and liver of the young mouse, and reverse similar signs in the old mouse 49,50 . Likewise, the detrimental cognitive effect of injecting old blood in mice is counteracted when anti-VCAM1 antibodies are concomitantly injected 40 . As such, the blood currently remains one of the most promising tissues to detect aging-related proteins.
Heme metabolism and iron levels were previously hypothesized to play a role in human aging 11 , and in the present study we identify the same pathway using new methods and additional data. It is of interest that both the heme pathway and the proteins we uncovered using MR are strongly linked to vascular and endothelial damage 51 . Across the genome, we also found an enrichment for brain tissues and pathways related to neuronal integrity. Given that endothelial cells are central to both the cardiovascular system and the blood-brain barrier 52,53 , and that endothelial function declines with age 54 , progressive endothelial dysfunction may manifest itself as an age-related disease. Indeed, recent findings suggest that the detrimental effects of the APOE ε4 allele-the largest genetic determinant of human aging-are mediated by an accelerated breakdown of the blood-brain barrier, independent of amyloid-β and tau accumulation 55 . We therefore speculate that molecules involved in maintaining or repairing endothelial integrity may be key to avoiding both age-related cardiovascular injury and neurodegeneration, and recommend further research into this area.
Our study demonstrates that GIP analysis of genetically correlated GWASs can increase power to detect shared genetic architecture and can identify genetic loci that would have been missed by any individual GWAS. For example, regions near AFF3, MAML3, USP28/HTR3B and MIR6074 reached genome-wide significance only in the combined analysis, and validation of MAML3 in an external survival cohort suggests that these findings may hold across populations and aging measurements. In addition, USP28 within the USP28/HTR3B locus resembles CDKN2A within the well-known CDKN2B/-AS1 locus, because both have been implicated in cellular proliferation 56,57 and senescence 58,59 . In fact, short hairpin RNA knockdown of USP28 in vitro results in decreased CDKN2A expression 58 . However, additional validation and functional investigation of aging-GIP1 loci are warranted.
A secondary advantage of the GIP method is the high stability of the resultant aging-GIP1 GWAS, which appears to be largely robust to the selection of component traits. Nevertheless, we note that our analysis focused on chronological and disease measures of aging and did not include the cellular measures of aging due to their modest genetic correlations and lack of power. Whether aging-GIP1 loci and pathways accurately reflect the aging process overall, or only capture specific domains of aging, is unknown. Inclusion of future aging-related GWASs performed on Biobank-scale samples should provide insight into how well aging-GIP1 captures human aging. An important limitation of the aging-GIP1 GWASs, and GWASs of aging-related traits more generally, is the potential confounding of aging genetics with social stratification. Socioeconomic factors exhibiting geographical clustering can induce gene-environment correlations that could inflate trait heritability and may bias results, especially for UK Biobank studies 25 . In our study, we found that genetic correlations of the six aging-related GWASs with socioeconomic deprivation were moderately high (up to 61% for self-rated health). This confounding can be inadvertently propagated and even amplified in the shared genetic component captured by aging-GIP1. Indeed, when explicitly removing genetic correlations with two UK Biobank GWASs of socioeconomic factors, aging-GIP1 heritability is halved. After adjustment, genetic correlations with mental and physical health traits are also somewhat reduced, although the overall patterns remain stable. Similarly, effect sizes of blood protein levels and aging-GIP1 loci-including loci that have been confidently linked to aging-related traits previously-were only slightly reduced after adjustment.
The shift of these effect sizes after adjustment was almost equal for all SNPs and MR signals. Together with the highly decreased heritability and stable genetic correlations, this suggests that the adjustment removed the unspecific polygenic background induced by UK Biobank microstructure and social confounding. We therefore consider this procedure to be a suitable approach for confounder correction, although it may be somewhat conservative for also removing the genetic background of the adjusting traits, which may genuinely be associated with aging. An alternative solution to the confounding problem could be the inclusion of a more ancestrally diverse set of aging-related GWASs from a larger variety of populations.
Despite these limitations, modeling the shared genetic component of human aging proxies has allowed us to downweigh nonbiological features, and propose pathways and proteins that may causally influence the human aging process. We share the full aging-GIP1 summary statistics without restrictions 60 to encourage further MR analysis using other biomarkers, and accelerate the discovery of drug targets able to prolong mental and physical wellbeing throughout life.

Methods
Data sources. We searched PubMed and Google Scholar in April 2020 for GWASs of aging measures, including only studies for which we could obtain autosomal genome-wide summary statistics measured in at least 10,000 European-ancestry individuals. If multiple studies were performed on similar traits, we kept the study with the largest sample size. GWASs meeting inclusion criteria included extreme longevity 17 (survival past 90th percentile), father and mother lifespan 4 , healthspan 61 , self-reported health 16 , frailty index 8 , epigenetic age acceleration 19 , telomere length 5 , mosaic loss of the Y chromosome 7 and perceived age 18 . As the self-reported health GWAS used the first release of UK Biobank data (n = ~150,000), we looked up the same phenotype in the Neale Lab GWAS collection (n = ~500,000), which had a larger sample size but was otherwise measured identically 62 . The original derivation of each set of summary statistics is briefly described in the Supplementary Note.
For each set of summary statistics, we discarded poorly imputed (INFO < 40%), rare (minor allele frequency (MAF) < 0.5%) or poorly measured SNPs (n individuals <1% of total). The remaining SNPs were aligned to genome build GRCh37 and harmonized to match UK Biobank SNP IDs (discarding any duplicates). We then estimated the phenotypic variance of the trait (residuals) from a selection of independent SNPs provided by the MultiABEL R package 13 where Var(Y) is the phenotypic variance, p and q are the major and minor allele frequencies, n is the total sample size (cases + controls, if applicable) and Var β is the variance of the effect size estimate. If n differed by SNP, we calculated the phenotypic variance separately for quintiles of n and took the mean estimate. SNP statistics were then standardized by dividing effect sizes and s.e. by √ Var(Y).  20 and the aging-related GWAS and aging-GIP summary statistics used in our study, where the longer computation time of the HDL software would have been impractical. The GWAS-MAP platform contains summary statistics for 1,329,912 complex traits and gene expression levels, and 3,642 binary traits, derived from UK Biobank 62,65 and large European-ancestry consortia (for example, MAGIC, CARDIoGRAM, SSGAC and GIANT; for an up-to-date list of phenotypes, see https://phelige.com) 66 . We included GWASs with ≥1 million SNPs, measured in ≥10,000 individuals (if continuous) or ≥2,000 cases and controls (if binary). We further excluded the healthspan and self-rated health GWASs from the GWAS-MAP platform to avoid duplication, after which 728 traits remained. Out of specific interest 67 , we also calculated the genetic correlation between aging-GIPs and a case-control GWAS of COVID-19 hospitalization, with cases defined as laboratory-confirmed COVID-19 patients experiencing a severe outcome and controls defined as the rest of the population (A2_ALL excluding 23andMe, release 4) 68 . GWASs with counterintuitive effect directions were reversed (for example, satisfaction traits are coded from high to low in UK Biobank), and each set of GWAS summary statistics was filtered before analysis by excluding SNPs in the MHC, SNPs with MAF < 1% and, if measured, SNPs with INFO < 90%. A full list of 729 GWASs and references can be found in Supplementary Data 3. P values were adjusted for multiple testing using Bonferroni's correction (729 traits and 6 GIPs).

Estimation of genetic and nongenetic correlations.
The same software was used to calculate pairwise correlations between GWAS-MAP statistics, which allowed traits to be clustered based on the magnitude of their genetic similarity. For computational tractability, we included only GWASs that showed large and significant effects on aging-GIP1 (|r g | > 0.25; P adj < 0.05). Clusters were identified hierarchically by maximizing the Bayesian information criterion using the mclust R package 69 v.5.4.1, up to a maximum of 100 clusters.

Identification of shared and unique genetic correlations.
For each selected GWAS-MAP phenotype, genetic correlations with aging-related traits were meta-analyzed using fixed-effect, inverse-variance weighting, with heterogeneity quantified by Cochran's Q and I 2 statistics, implemented in the meta R package 70 v.4.15-1. GWAS-MAP phenotypes that were significantly correlated with all six aging-related GWAS at FDR of 5% and lacking substantial heterogeneity (P het > 0.05 and I 2 < 50%) were considered to be shared. For GWAS-MAP phenotypes with evidence of heterogeneity, we performed a leave-one-out sensitivity analysis to assess which aging-related trait(s) contributed most to this heterogeneity. If heterogeneity could be completely removed by excluding a single GWAS (I 2 = 0%), and exclusion of any other GWAS did not substantially reduce heterogeneity (I 2 ≥ 50%), the aging-related trait outlier was considered to have a unique genetic correlation with the GWAS-MAP phenotype.
Genetically independent phenotype analysis of aging GWAS. Principal component loadings for six aging-GIPs were estimated from the genetic covariance matrix of the six chronological and holistic aging GWASs, analogous to a principal component analysis of phenotypic correlations. Specifically, eigenvectors from the genetic covariance matrix were transformed into loadings by dividing them by the square root of the phenotypic variance of the GIP. This phenotypic GIP variance was calculated as follows: where a i is the eigenvector of the ith aging-GIP and Σ ph is the phenotypic variancecovariance matrix of the core aging traits. As GWASs were standardized, Σ ph is equivalent to the phenotypic correlation matrix, with off-diagonal correlations estimated from the correlation observed between independent null z statistics (described above). The six chronological and holistic aging GWAS statistics were then combined on an SNP-by-SNP basis using the principal component loadings to construct genome-wide summary statistics for the six aging-GIPs. GIP effect estimates were calculated by summing effect estimates from the individual aging-related trait GWAS, each multiplied by their corresponding loading, and standardized using the expected GIP variance (equation (2)). The s.e.s of the GIP effect estimate were calculated by performing the equivalent calculation using variance arithmetic, also taking into account the phenotypic covariance between GWASs to adjust for sample overlap: where s.e.(B) is the vector of SNP s.e.s of the core aging trait effect estimates. Effective sample sizes were then estimated based on the median z statistic and allele frequencies, that is, solving equation (1) for n. Full technical details of the GIP method are described in Supplementary Note. Aging-GIP summary statistics were calculated for the 7,324,133 SNPs shared between quality-controlled GWAS statistics, of which 5,353,660 were common (MAF ≥ 5%) and 1,970,474 were rare (MAF < 5%). Finally, s.e.s of each aging-GIP GWAS were adjusted to account for the LDSC regression intercept, which ranged from 0.99 (GIP1) to 1.03 (GIP4).
Adjusting for genetic correlations with socioeconomic factors. UK Biobank GWAS summary statistics for household income and Townsend's deprivation index were downloaded from the Neale Lab GWAS collection 62 and subjected to the same SNP filters used for aging-related trait quality control. Genetic correlations between GWASs of aging-GIPs and the two socioeconomic GWASs were calculated using the HDL R package 21 v.1.3.4. As before, phenotypic correlations were estimated from the correlation between independent null z statistics.
For each aging-GIP, adjustment loadings were then calculated as: where β refers to β = Cov(X) −1 × Cov(GIP, X). Here, Cov(X) is the genetic covariance matrix of UK Biobank GWAS of household income and Townsend's deprivation index, and Cov (GIP, X) is the genetic covariance matrix between the aging-GIP and the two GWAS. These loadings were used as weights in the linear combination framework to calculate new aging-GIP summary statistics, which-by definition-were uncorrelated to the UK Biobank GWAS of household income and Townsend's deprivation index. Last, s.e.s were adjusted to account for the LDSC regression intercept, which ranged from 0.99 (GIP3) to 1.02 (GIP1).
Characterization of genome-wide significant aging-GIP1 loci. Genome-wide significant loci were defined as 500-kb regions centered on a lead genome-wide significant SNP (P < 8.3 × 10 −9 ) in linkage equilibrium (r 2 < 0.1) with other lead locus SNPs. The LD between SNPs was calculated using a random 10,000 subset of unrelated, genomically British individuals from UK Biobank 71 . The susieR package 26 v.0.11.8 was used to perform genetic fine-mapping of the association signals within each aging-GIP1 locus, allowing for up to ten causal variants. We report posterior inclusion probabilities for SNPs within 95% credible sets. Positional annotations of credible set SNPs were retrieved using ANNOVAR 72 v.2019Oct24.
SNPs within each 500-kb locus were looked up in the six aging-related trait GWASs used to construct the aging-GIPs. Any region containing an SNP associated with an aging-related trait at genome-wide significance (P < 5 × 10 −8 ) was tested for colocalization between the trait and the aging-GIP1 signal. For this, we used the coloc R package 32 v.4.0-4 with its default parameters, denoting a posterior colocalization probability (PP) of 80% as evidence of a shared GWAS signal.
Association of loci with Finnish and Japanese survival. Loci were considered to be previously replicated if they had been associated at genome-wide significance with one of the core aging traits and also had evidence of an effect in an independent cohort on the same trait (P < 0.05). We attempted to find additional evidence for an effect on aging for the aging-GIP1 loci, which had not been previously replicated.
Effects were first looked up in a GWAS of survival of FinnGen study participants 73 . The FinnGen study associated SNPs across the genome with the survival of 218,396 Finnish-ancestry individuals (203,244 censored, 15,152 deceased) using Genetic Analysis of Time-to-Event phenotypes (GATE) v.0.40. SNP effects were log(hazard ratios), calculated from a mixed-effect frailty model that adjusted for sex, genotyping batch, birth year and the first ten genomic principal components as fixed effects, and cryptic relatedness using the genetic relatedness matrix as random effects.
Analogously, the same SNPs (if polymorphic) were regressed against the survival of 135,983, unrelated, Japanese-ancestry individuals (97,365 censored, 30,976 deceased) from BioBank Japan 74,75 . In this analysis, a fixed-effect Cox's proportional hazards model was fitted using the survival R package v.2.41: where h(t) is the hazard at time t, given that the subject is alive at time t and h 0 (t) is the baseline hazard at time t; X 1 , X 2 , …, X n are the vectors of covariates with fixed effects β 1 , β 2 , …, β n ; and γ is the effect of the vector of SNP dosages G. Covariates were sex, disease status and the first 20 principal components, where disease status refers to 1 of 47 common diseases in Japan used to recruit the individuals. Each SNP was tested in a separate model. The SNP effects from both studies were converted from log(hazard ratios) to approximate years of life by inverting the sign and multiplying the effect estimate and s.e.s by ten 3 . For each SNP, a combined effect was calculated by meta-analyzing the cohort-specific effects in a fixed-effect framework (weighted using inverse variance), implemented in the meta R package 70 v.4.15-1. One-sided P values were adjusted for multiple testing of 23 loci using Bonferroni's correction. The collective effect of the 21 remaining loci was calculated using a random-effect framework, to allow for heterogeneity in effect size estimates.
Aging-GIP1 leave-one-out sensitivity analyses. Leave-one-out sensitivity analyses were performed for aging-GIP1, where, one at a time, a core aging trait was excluded and aging-GIP1 loadings and summary statistics were recalculated using the remaining five traits. Genetic and nongenetic correlations were calculated between the original aging-GIP1 and each leave-one-out GIP using HDL inference and null SNP z statistics, as described above.
To test for heterogeneity in genome-wide significant aging-GIP1 loci, we estimated the difference between the lead SNP effect in aging-GIP1 and the effect in the leave-one-out GIP1 GWAS, taking into account null correlations between traits. The s.e. of the difference in effects was calculated as follows: s.e. β1 − β2 = s.e. β1 2 + s.e. β2 2 − 2r × s.e. β1 × s.e. β2 (5) where s.e. (β1) and s.e. (β2) are the s.e.s of the SNP for aging-GIP1 and the leaveone-out GIP1, respectively, and r is the nongenetic correlation between GWASs. Statistical significance of the difference was assessed using Wald's test and adjusted for multiple testing of 27 loci using Bonferroni's correction.
Lookup of known SNP associations. Lead SNP and close proxies (r 2 EUR ≥ 0.8) of the aging-GIP1 loci were looked up in PhenoScanner 27 and the GWAS catalog 28 (accessed 3 December 2020), keeping only the traits with genomewide significance (P < 5 × 10 −8 ). Triallelic SNPs and associations with treatments or medications were discarded, before converting associations with the lack of a phenotype into the phenotype itself by inverting the sign (for example, 'Qualifications: none' to 'Qualifications'). We then further grouped the traits based on similarities in trait names, keeping the strongest association in the group. This grouping was done by partial matching of trait names-verified manually-and keeping the shortest name. For example, 'Melanoma' , 'Malignant melanoma' and 'Malignant melanoma of skin' were grouped and renamed to 'Melanoma' .
Tissue enrichment. Stratified LDSV regression v.1.0.0 was used to stratify aging-GIP1 SNPs into categories and test whether the proportion of SNP heritability in a category exceeded that expected from the proportion of SNPs in the category 64 . We kept only HapMap3 SNPs, excluding the MHC region and SNPs with MAF < 0.05, and used the 1000 Genomes Phase 3 LDSC reference as weights. Categories tested included the 10 groups summarizing 220 cell-type specific annotations from Finucane et al. 64 , adjusting for the baseline model (v.1.2).
Gene and pathway enrichment. PASCAL 29 was used to aggregate aging-GIP1 SNP-level P values (with and without adjustment for socioeconomic correlations) into gene scores and test these scores for enrichment against predefined gene sets. Gene sets were Hallmark (C1) and Gene Ontology Biological Process (C5.BP) sets from v.7.2 of the Molecular Signatures Database 76 .
Aging-GIP1 summary statistics were first aligned to the 1000 Genomes SNP build (matching the PASCAL LD reference), before being tested with default PASCAL parameters, which includes discarding SNPs with MAF < 5% and SNPs in the MHC region. Gene results passing a 5% FDR threshold were considered to be significant. For each pathway in the C1 and C5.BP datasets, PASCAL calculated two measures of significance based on χ 2 and permutation statistics. We separately adjusted C1 and C5.BP for multiple testing and considered a pathway with both χ 2 and permutation statistics passing a 5% FDR threshold to be significant.
Significant C5.BP pathways (with and without adjustment for socioeconomic correlations) were clustered based on their Jaccard's similarity coefficient (that is, the size of the intersection of genes divided by the size of the union of genes). The mclust R package 69 v.5.4.1 was used to maximize the Bayesian information criterion of Jaccard's similarity matrix to identify the optimal number of clusters (up to a maximum of 100). We used the number of clusters selected by most mclust models to group the pathways.
MR of blood protein levels. Genetic instruments for blood protein levels (pQTLs) were retrieved from Zheng et al. 30 . We included all tier 1 instruments: cis-and trans-pQTLs shown to influence ≤5 proteins (specificity) with no evidence of heterogeneity in effect sizes between multiple protein expression studies (consistency). A total of 857 proteins had nonpalindromic SNP instruments (898 total pQTLs) present in the aging-GIP1 summary statistics. Two-sample MR of blood protein levels as exposures and aging-GIP1 as outcome was performed using the TwoSampleMR R package 31 v.0.5.5. If multiple pQTL instruments were available for a protein, heterogeneity and MR-Egger sensitivity tests were also performed. This analysis was repeated with the six standardized aging-GIP1 component traits as outcome, as well as the unstandardized, combined parental lifespan GWAS from Timmers et al. 4 to provide an intuitive measure of the effect.
The MR effects and s.e.s from the latter were multiplied by 10 to convert them from units of negative log(hazard ratio) to approximate years of life 3 .
Aging-GIP1 MR results passing a 5% FDR threshold and sensitivity tests (P Het > 0.05 and P Egger > 0.05, if applicable) were taken forward for follow-up colocalization tests to rule out LD linkage. For proteins for which we had access to genome-wide summary statistics (for example, LPA, β 2 M and VCAM1), we used the coloc R package 32 v.4.0-4 to perform colocalization analysis using the default parameters and denoted PP ≥ 80% as evidence of a shared signal. For instruments without summary statistics, we performed an LD check as described in Zheng et al. 30 , which involved checking whether any of the 30 strongest aging-GIP1 SNPs in a 1-Mb region centered on each pQTL were in high LD (r 2 EUR ≥ 0.8) with that pQTL.
Finally, proteins passing sensitivity and colocalization tests were subjected to a reverse-causality test using MR-Steiger 77 and, if genome-wide summary statistics were available, a bidirectional MR analysis, both implemented in TwoSampleMR. For the bidirectional MR, we used up to 27 genome-wide significant lead SNPs from aging-GIP1 (shared between GWASs and replacing missing or palindromic SNPs with the next most significant SNP) as instruments and the protein expression statistics as outcome. Proteins significant for the MR-Steiger test (P < 0.05) and showing no evidence of reverse causality in the bidirectional MR (P > 0.05), if applicable, were considered to have robust causal effects on aging-GIP1.
Ethical oversight. BioBank Japan participants provided written informed consent and survival study protocols were approved by BioBank Japan Project ethical review boards from the Institute of Medical Sciences, University of Tokyo and the RIKEN Center for Integrative Medical Sciences. All other data were publicly available and approved by ethical committees as described in their respective publications. See Supplementary Note for FinnGen ethical approval.

Statistics and reproducibility.
The statistical methods used to analyze the data are described fully in Methods, with basic data processing done using R v.3.6.0 unless specified otherwise. We used the independent FinnGen and BioBank Japan cohorts to successfully replicate 2 of the 23 new aging-related loci. Predetermination of study sample size, randomization of experiments and blinding of investigators to experiments were not applicable for our type of study.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.