Genome-wide association studies identify 137 genetic loci for DNA methylation biomarkers of aging.

BACKGROUND: Biological aging estimators derived from DNA methylation data are heritable and correlate with morbidity and mortality. Consequently, identification of genetic and environmental contributors to the variation in these measures in populations has become a major goal in the field.
RESULTS: Leveraging DNA methylation and SNP data from more than 40,000 individuals, we identify 137 genome-wide significant loci, of which 113 are novel, from genome-wide association study (GWAS) meta-analyses of four epigenetic clocks and epigenetic surrogate markers for granulocyte proportions and plasminogen activator inhibitor 1 levels, respectively. We find evidence for shared genetic loci associated with the Horvath clock and expression of transcripts encoding genes linked to lipid metabolism and immune function. Notably, these loci are independent of those reported to regulate DNA methylation levels at constituent clock CpGs. A polygenic score for GrimAge acceleration showed strong associations with adiposity-related traits, educational attainment, parental longevity, and C-reactive protein levels.
CONCLUSION: This study illuminates the genetic architecture underlying epigenetic aging and its shared genetic contributions with lifestyle factors and longevity.

Background

Aging is associated with an increased risk of physical, cognitive, and degenerative disorders [1]. While the rate of chronological aging is constant between individuals, there are inter-individual differences in the risk of age-associated morbidities. Biological aging is influenced by both environmental and genetic factors [2]. Multiple measures of biological age exist, several of which have drawn information from DNA methylation (DNAm) across the genome. DNAm is a common epigenetic modification typically characterized by the addition of a methyl group to a cytosine-guanine dinucleotide (CpG). DNAm levels can be influenced by both genetic and environmental factors, and in recent years, DNAm signatures have become established correlates of multiple health-related outcomes [3-5]. Such signatures include “epigenetic clocks”, accurate markers of aging which associate with several health outcomes [6, 7]. Epigenetic clocks use weighted linear combinations of CpGs to predict an individual’s chronological age and have common single-nucleotide polymorphism (SNP)-based heritability estimates ranging from 0.15 to 0.19 [8, 9]. Individuals with epigenetic clock estimates greater than their chronological age display “age acceleration” and have been shown to be at a greater risk of all-cause mortality and multiple adverse health outcomes [10]. Consequently, identification of genetic and environmental contributors to the variation in these measures in populations has become a major goal in the field [11].

The first generation of epigenetic aging clocks used penalized regression models to predict chronological age on the basis of DNA methylation data, e.g., the widely used clocks from Hannum (2013) and Horvath (2013) apply to blood and 51 human tissues/cell types, respectively [12-14]. A derivative of the Horvath clock, intrinsic epigenetic age acceleration (IEAA) has since been developed, conditioning out (i.e., removing) estimates of blood cell composition. An increasing literature supports the view that IEAA relates to properties of hematopoietic stem cells [2, 8, 15]. The second generation of epigenetic clocks move beyond estimating chronological age by incorporating information on morbidity and mortality risk (e.g., smoking, plasma protein levels, white blood cell counts), and chronological age. Two such predictors, termed PhenoAge (a DNAm predictor trained on a measure that itself was trained on mortality, using 42 clinical measures and age as input features) and GrimAge (trained on mortality, including a DNAm measure of smoking as a constituent part), outperform both Hannum and Horvath clocks in predicting mortality and are associated with various measures of morbidity and lifestyle factors [16, 17]. DNAm GrimAge outperforms PhenoAge and the first generation of epigenetic clocks when it comes to predicting time to death [8, 18, 19].

While nothing is known about the genetics of the second generation of epigenetic clocks, 13 genetic loci have been associated with the first generation of epigenetic clocks. A study of nearly ten thousand individuals revealed a regulatory relationship between human telomerase (hTERT) and epigenetic age acceleration [8]. More recently, a larger genome-wide association study (GWAS; n = 13,493) revealed that metabolic and immune pathways share genetic underpinnings with epigenetic clocks [9].

Here, we greatly expand on these studies across several dimensions. First, we analyze a large, multi-ethnic dataset comprised of over 41,000 individuals from 29 European ancestry studies, seven African American studies, and one Hispanic ancestry study. Second, we characterize for the first time the genetic architecture of the second-generation epigenetic clocks, GrimAge and PhenoAge. All of these clocks have been trained on European ancestry populations. Third, we also conduct GWAS of two important DNAm-based surrogate markers: DNAm plasminogen activator inhibitor-1 (PAI1) levels and granulocyte proportion, respectively. Although not considered in risk prediction scores such as the Framingham Heart score, DNAm PAI1 was chosen because it exhibited stronger associations with cardiometabolic disease than the epigenetic clocks [16]. The DNAm-based estimate of granulocyte proportions was chosen because it exhibited significant associations with several epigenetic clocks (including GrimAge and PhenoAge) and with health outcomes such as Parkinson’s disease [16, 17, 20]. The unprecedented sample size of the current study allowed us to develop polygenic risk scores for these six epigenetic biomarkers.

We report 137 independent loci, including 113 novel loci (i.e., not previously identified in previous GWAS meta-analyses of epigenetic age estimators [8, 9]), and examine the genetic and causal relationships between epigenetic aging, lifestyle behaviors, health outcomes, and longevity.

Results

To identify genetic variants associated with six methylation-based biomarkers, genome-wide association studies of 34,710 European ancestry and 6195 African American individuals were performed (Additional file 1 and Additional file 2: Tables S1-S2). A fixed effects meta-analysis was performed to combine the summary statistics within each ancestry group (summary statistics available at https://datashare.is.ed.ac.uk/handle/10283/3645). Genomic inflation factors ranged between 1.01 and 1.06 (Additional file 3: Figure S1-S6) for the European-only meta-analyses, indicating appropriate adjustment for population stratification, and from 1.11 and 1.21 for the meta-analyses comprising African American participants (Additional file 3: Figures S7-S12; Table 1). Inflation was present and consistent across all allele frequencies in the African American analyses; there was much greater variability in the effect sizes in the African American analyses (Additional file 4). Phenotypic correlations were examined in Generation Scotland, the largest participating cohort in the study. Correlations ranged from 0 (IEAA and granulocyte proportions) to 0.48 (PhenoAge acceleration and Hannum age acceleration (Additional file 2: Table S3). We examined the relationship between means and standard deviations of predicted age versus means and standard deviations of chronological age for each cohort, separated by ancestry group, observing weaker mean correlations in the African American cohorts. There was little difference in the relationship between the standard deviations of age acceleration and chronological age by ancestry group (Additional file 3: Figure S13). Heterogeneity between studies may decrease power to detect genetic associations. We found little evidence of systematic between-study heterogeneity in both the European ancestry and African American meta-analyses, as determined by M-statistic outlier analysis, meta-regressions against cohort characteristics, and analysis of heterogeneity I2 statistics [21] (Additional file 4).

Table 1

Summary of key findings of post GWA analyses. African American meta-analyses based on 6195 participants from 7 cohorts; European ancestry meta-analyses based on 34,710 participants from 28 cohorts

The key findings of the post GWA analyses are summarized in Table 1 along with a summary of the input features of each clock in Table 2. The latter highlights the value in discriminating novel associations from those that are likely driven by the construction of the clocks.

Table 2

Summary of input features of first-generation (Hannum and Horvath) and second-generation (DNAm PhenoAge, GrimAge) epigenetic clocks

European ancestry GWAS meta-analysis: 56 independently associated loci

We identified 56 conditionally independent associations (P < 5 × 10−8) across the six epigenetic biomarkers in European ancestry populations using a stepwise model (Additional file 3: Figures S14-S19; Table 1; Additional file 2: Table S4) [22, 23]. We replicated 10/10 loci associated with IEAA and 1/1 locus associated with a cell-adjusted Hannum-based measure of epigenetic age acceleration identified in an earlier GWAS (P < 0.05/11 = 0.0045) [9]. All but three loci (associated with IEAA) were replicated at the genome-wide significant level. To validate the associations with DNAm-derived granulocyte counts, we compared our results to a previous GWAS of FACS granulocyte counts which identified 155 independent loci [24]. In the current meta-analysis, we replicated 13/129 present loci (P < 0.05/129 = 3.88 × 10−4; Additional file 2: Table S5; Additional file 3: Figures S20-S21), two of which replicated at the genome-wide significant level. Effect sizes at the 129 loci were strongly correlated between studies (r = 0.85). Conversely, we failed to replicate four genome-wide significant lead SNPs from a previous GWAS of measured PAI1 levels (P ≥ 0.327) [25]. There was no clear concordance of effect sizes at these loci in the current study (r = -0.59; P = 0.41)

To examine whether genetic variation across the six epigenetic biomarkers was shared, we performed genetic colocalization analyses of the 56 loci [26]. There was evidence for colocalization at 30 loci between epigenetic biomarkers (posterior probability (PP) > 0.8; Additional file 2: Table S6). IEAA was associated with the greatest number of independent loci (n = 24), whereas granulocyte proportion was associated with the fewest (n = 2).

African American GWAS meta-analysis: 81 independently associated loci

We identified 81 conditionally independent associations for the six epigenetic biomarkers in the African American analyses (Additional file 3: Figures S22-S27; Table 1; Additional file 2: Table S4). The number of associated loci per epigenetic biomarker ranged from 9 (IEAA) to 27 (granulocyte proportion).

Trans-ethnic meta-analyses identify 69 loci

To determine if any loci were shared across the European and African American populations, a trans-ethnic meta-analysis was carried out for each of the six epigenetic biomarkers using MR-MEGA [27]. Sixty-nine risk loci were identified across the six predictors that were common to all ancestries, ranging from 6 (GrimAge acceleration) to 23 (IEAA). Ten loci were significant in the African American analyses, 33 were significant in the European analyses, and five were significant in both. This left 21 novel loci from the trans-ethnic meta-analyses (Additional file 2: Table S7). Among the allele frequencies of the lead SNPs for these loci, 11/21 differed by > 10% between European and African American populations (Additional file 2: Table S8).

We compared effect sizes of the lead SNPs from these loci in a Hispanic-American ancestry subset of the MESA cohort (n = 287). Correlations between the respective effect sizes ranged from very weak (r = 0.16 for 10 granulocyte proportion SNPs) to near unit (r = 0.92 for 10 Hannum age acceleration SNPs; Additional file 2: Table S9 and Additional file 3: Figure S28).

Gene-based GWAS identifies 364 significant genes

Gene-based GWASs carried out using MAGMA identified between two and 46 genes (111 unique genes in total) associated with the six epigenetic biomarkers in the European ancestry data (Additional file 2: Table S10) [28]. In the African American data, between nine and 209 genes (264 unique genes in total) were associated with the epigenetic biomarkers (Additional file 2: Table S10). Across all epigenetic biomarkers and ancestries, there were 364 unique gene-based associations.

Independently associated loci are associated with DNA methylation levels

One obvious genetic effect that may influence our GWA findings is the overlap with cis methylation quantitative trait loci (mQTLs) for epigenetic clock DNAm sites. To explore whether any of the 56 loci from the European GWAS shared genetic variation influencing epigenetic clock DNAm sites, colocalization analyses were conducted using GoDMC summary statistics [29] (“Methods”). We found strong evidence (PP > 0.8) that 1/4 loci (25%) for GrimAge acceleration, 3/12 loci (25%) for PhenoAge acceleration, 11/24 loci (46%) for IEAA, and 5/9 loci (56%) for Hannum age acceleration had shared genetic variation influencing epigenetic clock DNAm sites (Additional file 2: Tables S11-S12). Next, we used genetic colocalization to assess whether GWAS loci for aging biomarkers were associated with methylation levels at established DNAm sites for BMI [30] and smoking [31]; 29/56 loci (52%) colocalized with genetic variation influencing smoking-associated DNAm sites and 1/56 loci (1.8%) was colocalized with genetic variation influencing a BMI-associated DNAm site. Specifically, GrimAge acceleration (75%) and Hannum age acceleration (78%) loci showed a large overlap of genetic variation influencing smoking DNAm sites.

Utilizing results from a published GWAS of IEAA in brain tissue [32], we tested whether genetic variation influencing IEAA in blood and brain was shared for the 24 blood-related IEAA loci (cross-tissue plot for lead SNPs shown in Additional file 3: Figure S29). Colocalization analysis showed that there was no strong evidence (PP > 0.80) for a single SNP being associated with both traits. However, the true extent of sharing is difficult to estimate because the sample size of the brain study (n = 1796) is much smaller than our blood-based study, limiting power to detect shared loci. Previous simulations using a sample size of 2000 individuals have indicated that the shared variant must explain close to 2% of the variance of a biomarker to attain a posterior probability > 0.8 for shared genetic effects [26]. Nevertheless, we observed suggestive evidence for colocalization (PP = 0.53; LocusZoom plot in Additional file 3: Figure S30) for a locus mapping to DSCR6 on chromosome 21. This locus also shares genetic variation with an IEAA clock CpG, cg13450409 (PP = 0.99; Additional file 2: Table S12).

SNP- and gene-based enrichment within published GWAS

To determine whether any of the 56 lead SNPs in the European ancestry meta-analyses for the six epigenetic biomarkers showed evidence for pleiotropic associations, a lookup of published GWAS significant associations (P < 5 × 10−8) was carried out (Additional file 2: Table S13). Four of the IEAA-associated SNPs (rs2736100 in TERT, rs2275558 in PBX1, rs144317085 in TET2, and rs2492286 in RPN1) were associated with 16 unique traits including multiple cancers (e.g., lung cancer, glioma) [33, 34] and blood cell counts (e.g., platelet count, eosinophil count, red blood cell count) [24, 35]. Whereas there was considerable overlap with cell-related traits (e.g., white/red blood cell, platelet, and eosinophil counts, and mean corpuscular volume) [24, 35, 36], there were no associations with non-cancer-related disease or lifestyle measures (Additional file 2: Table S13). A gene-based test of enrichment among traits within the GWAS catalog output (Additional file 2: Table S14) showed genes associated with IEAA, GrimAge, and granulocyte proportion were enriched among those associated with white blood cell counts [24, 35]. Several genes associated with granulocyte proportions were also enriched among those associated with inflammatory traits (e.g., inflammatory bowel disease, rheumatoid arthritis, asthma) [37-42]. IEAA-associated genes were also significantly enriched among those identified in a previous GWAS of IEAA [8, 9].

Colocalization to identify GWAS loci that might regulate expression levels

Colocalization analyses were conducted to investigate whether any of the 56 loci from the European ancestry GWASs showed evidence of regulating gene expression levels (“Methods”). There was strong evidence (PP > 0.8) that eight loci had shared genetic effects with expression quantitative trait loci (eQTLs; Table 1; Additional file 2: Table S15). Of these, one was associated with GrimAge acceleration and seven were associated with IEAA. The locus associated with GrimAge acceleration was linked to the expression of C6orf183 whereas IEAA-associated loci were linked to the expression of 11 transcripts including genes related to lipid transport and immune function (e.g., ATP8B4, CD46, TRIM59). Colocalization plots are presented in Additional file 5. Notably, four loci (all associated with IEAA) were independent of variants associated with DNAm levels at constituent clock CpGs (Additional file 2: Table S12; Table S15).

Functional enrichment analysis

To gain an understanding of the regulatory properties of the variants that underlie the six epigenetic biomarkers, we performed functional enrichment analyses across various gene annotations and regulatory and cell-type specific elements on the summary statistics for each of the European ancestry GWAS results (see “Methods,” Additional file 2: Table S16) [43]. At an epigenetic biomarker-specific adjusted P value calculated from the effective number of annotations, significant enrichments were present for IEAA (n = 191) and Hannum age acceleration (n = 3). Associations with IEAA were enriched in DNaseI hypersensitive site (DHS) hotspots in several tissues, which might reflect that Horvath’s pan-tissue clock applies to all tissues. The strongest enrichment of associations with IEAA could be observed for mobilized CD34 primary cells (OR = 6.06, P = 6.1 × 10−12), which supports the view that IEAA reflects properties of hematopoietic stem cells [2, 8, 15].

Pathway enrichment analysis

In the African American analysis, genes associated with granulocyte proportions were enriched among 35 Gene Ontology (GO) terms (Bonferroni P < 0.05), the majority of which were immune-related (e.g., adaptive immune response, lymphocyte activation, regulation of immune system process, or skin development; Additional file 2: Table S17). By contrast, there was no significant enrichment among Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways or GO terms for the significantly associated genes in the European ancestry analysis.

Overlap with Mendelian disease genes

In order to characterize the potential overlap of Mendelian disease genes and associated pathways with our findings, a series of enrichment analyses were conducted (“Methods”) [44]. In the European ancestry analysis, enrichment of Mendelian disease genes was observed for IEAA and two gene sets (“disorders of platelet function” and “vascular skin abnormality”: bootstrapped P = 0.027 and 0.049, respectively). Enrichment analysis of PhenoAge loci revealed overrepresentation of methylation-related Mendelian disease genes (Additional file 2: Table S18; bootstrapped P = 0.04) with MTR (methyltetrahydrofolate-homocysteine S-methyltransferase) present in addition to TPMT. No significant over-representations of Mendelian disease gene sets were observed for any of the genes identified in the African American analysis.

Heritability and LD score regression

We quantified the proportion of variance in the six epigenetic biomarkers from the European ancestry meta-analyses that can be explained by our SNP sets using Linkage Disequilibrium (LD) Score regression [45]. The GWAS summary statistic SNP-based heritability ranged from 0.02 (SE = 0.02) for PAI1 levels, to 0.17 (SE = 0.02) for IEAA (Table 1; Additional file 3: Figure S31; Additional file 2: Table S19). We omitted DNAm PAI1 from the genetic correlation analysis due to its low heritability estimate. Several of the remaining five epigenetic biomarkers exhibited significant (P < 0.05) pairwise genetic correlation coefficients ranging from r = 0.28 (GrimAge acceleration and IEAA) to r = 0.66 (GrimAge acceleration and PhenoAge acceleration; Fig. 1; Additional file 2: Table S20).

Fig. 1

Genetic correlations and standard errors for the epigenetic biomarkers. IEAA (intrinsic epigenetic age acceleration)

For each epigenetic biomarker, large-scale genetic correlation analyses were conducted with 693 different traits. A selection of significant associations after Bonferroni correction for multiple testing (P < 0.05/693 = 7.22 × 10−5) is presented in Fig. 2.

Fig. 2

Significant genetic correlations between GrimAge acceleration (A) and PhenoAge acceleration (B) and a selection of GWAS traits. *This variable was originally coded with a high score representing lower health rating. We have multiplied the genetic correlation by − 1 for interpretability

PhenoAge acceleration had significant genetic correlations with 36 health-related traits including educational and cognitive traits (e.g., years of schooling, intelligence; rg = − 0.26 and − 0.30; P ≤ 3.29 × 10−5) [46-48], anthropometric traits (e.g., waist circumference, obesity, extreme BMI, hip circumference; rg = 0.22-0.31; P ≤ 3.61 × 10−5) [49, 50], (http://www.nealelab.is/uk-biobank/), adiposity (e.g., leg, arm, and trunk fat mass; rg = 0.22–0.23; P ≤ 3.46 × 10−6) (http://www.nealelab.is/uk-biobank/), and longevity (e.g., father’s age at death: rg = − 0.34; P = 9.66 × 10−6) (http://www.nealelab.is/uk-biobank/). GrimAge acceleration was genetically correlated with similar traits to PhenoAge (e.g., father’s age at death: rg = − 0.64; P = 6.2 × 10-13), along with smoking-related traits (e.g., current tobacco smoking: rg = 0.62; P = 1.5 × 10−15) (http://www.nealelab.is/uk-biobank/) and cancer-related traits (e.g., lung cancer: rg = 0.48; P = 8.3 × 10−6) [51]. The shared genetic contributions to GrimAge and smoking/mortality are expected given that GrimAge uses a DNAm-based estimator of smoking pack-years in its definition. There were no significant genetic correlations between Hannum age acceleration, IEAA, or granulocyte proportions and any of the traits tested after correction for multiple testing (Additional file 2: Table S21).

Polygenic risk score (PRS) profiling

To determine how well SNP-based genetic scores can approximate the six epigenetic biomarkers and investigate whether these genetic scores associate with health outcomes, a polygenic risk score analysis was conducted on the European ancestry data. Re-running the meta-analysis with an iterative leave-one-cohort-out process (and on the full summary statistics in a completely independent cohort—the Young Finns Study), the mean polygenic predictions explained between 0.21 and 2.37% of the epigenetic biomarkers (Table 1; Fig. 3; Additional file 2: Table S22). The maximum prediction for a single cohort was 4.21% for PAI1 levels in ARIES. Parsimonious predictors (built using SNPs with P < 5 × 10−8) performed well for IEAA, PAI1 levels, and PhenoAge acceleration, whereas predictors including more SNPs (P < 0.01–P < 1) tended to explain the most variance in GrimAge acceleration and granulocyte proportions.

Fig. 3

A Polygenic predictions for the six epigenetic biomarkers in LBC1921, LBC1936, SABRE, Born in Bradford, ARIES, FHS, and the Young Finns Study. IEAA (intrinsic epigenetic age acceleration), PAI1 (plasminogen activator inhibitor-1). B Associations between GrimAge polygenic risk score (P < 1) and UK Biobank GWAS traits

In order to investigate the association between the polygenic risk scores and health outcomes, we utilized a PRS Atlas to model associations with 581 heritable traits (nrange = 10,299 to 334,915) from the UK Biobank study [52]. The PRS inputs included the independent SNPs with P < 1 for GrimAge acceleration and granulocyte proportion and P < 5 × 10−8 for the other four epigenetic biomarkers, with thresholds based on the results from the leave-one-out predictions (Fig. 3). Using a false discovery rate (FDR)-corrected P value for each of the six epigenetic biomarkers, we found between 7 and 250 significant associations for GrimAge acceleration, granulocyte proportions, and Hannum age acceleration (PFDR < 0.05; Additional file 2: Table S23). The strongest associations were between the GrimAge acceleration PRS and the following traits: adiposity-related traits (e.g., body fat percentage: β = 0.02; PFDR = 7.3 × 10−39); education (e.g., college or university degree: β = − 0.06; PFDR = 2.6 × 10−43); and parental longevity (e.g., father’s age at death: β = − 0.02; PFDR = 5.7 × 10−16; mother’s age at death: β = − 0.02; PFDR = 1.6 × 10−11). Higher C-reactive protein was associated with a higher PRS for both granulocyte proportions and GrimAge acceleration (granulocyte proportions: β = 0.01; PFDR = 8.2 × 10−4; GrimAge acceleration: β = 0.02; PFDR = 2.1 × 10−29), and a lower score for Hannum age acceleration (β = − 0.006; PFDR = 0.02). A higher Hannum age acceleration PRS was also associated with an increased likelihood of taking insulin medication and lower total protein levels.

Mendelian randomization between age acceleration phenotypes and health and lifestyle outcomes

To investigate if the epigenetic measures were causally influenced by lifestyle factors and had a causal effect on aging and disease outcomes, we performed Mendelian randomization (MR) analyses on 150 traits for the European ancestry data (Additional file 2: Table S24). We found 12 inverse-variance weighted MR effects between the main exposures and epigenetic outcomes (GrimAge acceleration, PhenoAge acceleration, and PAI1 levels), after adjustment for multiple testing (Table 1; Additional file 2: Table S25). Of these, three remained significant (P < 0.05) across the other three MR methods. All of these consistent effects were with GrimAge as the outcome. Greater adiposity was associated with greater GrimAge acceleration: body mass index (BMI; BetaIVW = 0.76 years per standard deviation (SD) increase in BMI, P = 3.7 × 10−16); hip circumference (BetaIVW = 0.42 years per SD increase in hip circumference, P = 2.5 × 10−5); waist circumference (BetaIVW = 0.59 years per SD increase in waist circumference, P = 5.9 × 10-6; Fig. 4). Current tobacco smoking showed evidence for a causal effect on increased GrimAge in two of the MR methods (BetaIVW = 3.42 years for smokers, P = 9.0 × 10−6; Fig. 4), as anticipated given that it incorporates a DNAm-based estimator of smoking pack-years [8]. Past tobacco smoking showed evidence for an inverse causal effect (BetaIVW = − 1.09 years, P = 6.6 × 10−9), indicating that GrimAge acceleration is reduced upon smoking cessation. As a DNAm-proxy for leptin was included in the derivation of GrimAge, the smoking and adiposity findings may act as positive controls. There was evidence from three of the four MR methods to support a link between higher educational attainment (both years of schooling and college/university degree) and lower GrimAge acceleration (Fig. 4). For the secondary exposures, there was evidence across all methods for a causal effect of a greater body size at age 10 on higher GrimAge acceleration (BetaIVW = 0.70, P = 1.6 × 10−4; Additional file 2: Tables S24; S26). Consistent findings across all four MR methods provided evidence to support a causal effect of 13 cell count traits and DNA methylation-estimated granulocyte proportions, and between lower lymphocyte proportions and higher Hannum and GrimAge acceleration (Additional file 2: Table S27). There was evidence for heterogeneity in the causal effects for most of the cell types on epigenetic age measures, as well as years of schooling on GrimAge acceleration, although weaker evidence for directional pleiotropy was detected, based on the Egger intercept (Additional file 2: Table S28).

Fig. 4

Causal effects of UK Biobank GWAS traits on GrimAge acceleration. Effects correspond to increase/decrease in GrimAge acceleration per SD increase in waist circumference, hip circumference, and BMI; or per log odds increase for university/college education and current smoker status

The biomarker analyses with epigenetic measures as outcomes identified no consistent effects across all four MR methods (Additional file 2: Table S29). We found limited evidence to support any causal effects of the epigenetic measures (as exposures) on key disease and health outcomes, including longevity (Additional file 2: Tables S30-S31).

Discussion

Epigenetic biomarkers of aging and mortality have been extensively studied in relation to a plethora of health and disease outcomes. Here, we conducted a comprehensive suite of analyses in a meta-analysis sample of over 40,000 individuals, including the first GWA studies of two DNAm-based estimators of mortality risk (PhenoAge and GrimAge), as well as for DNAm-based proxies for granulocyte proportion and plasminogen activation inhibitor 1 protein levels. We identified 137 loci, of which 113 were novel, related to the six epigenetic biomarkers. Whereas previous studies have shown a general decline in longitudinal Hannum and Horvath age acceleration [53], there was no evidence of heterogeneity by cohort age for the lead independent loci. Although our comparison of genetic architectures across different ancestries was limited by sample size, particularly in the Hispanic ancestry lookup cohort of just 287 individuals, our European and African American trans-ethnic meta-analyses implicated many shared genetic loci. However, heterogeneity of effect sizes between European and African American ancestries was found for genome-wide significant loci which may reflect differential tagging of underlying causal variants. Alternatively, gene-environment interactions or poor prediction of epigenetic biomarkers in African Americans, possibly due to the construction of the clocks relying on European ancestry individuals, may explain heterogeneity in these effect sizes.

The IEAA clock is a derivative of Horvath’s pan-tissue predictor that regresses out DNA methylation-based estimates for naive CD8+ T cells, exhausted CD8+ T cells, plasmablasts, CD4+ T cells, natural killer cells, monocytes, and granulocytes. Although other cell sub-types may influence the findings, the GWAS results that are shared between IEAA and the other epigenetic clocks (Hannum clock, PhenoAge, GrimAge) are less likely to be influenced by differences in cell composition. However, there was still an enrichment in the GWAS results for IEAA and GrimAge with findings from previous white blood cell GWASs, such as basophils. Mendelian randomization analyses indicated putative causal effects of lymphocyte count on PhenoAge, GrimAge, and Hannum age acceleration but not on IEAA. Furthermore, whereas we observed enrichment of Mendelian disease genes for IEAA relating to platelet disorders, the Mendelian randomization analyses did not support evidence for a causal link between platelets and IEAA. IEAA and PhenoAge acceleration share the following genome-wide significant gene-based associations: TPMT, TERT, NHLRC1, KDM1B, EDARADD. IEAA and Hannum age acceleration share the following associations: TERT, TRIM59, KPNA4, RP11-432B6.3, IFT80, and TET2. TET2 is particularly interesting in light of its mechanistic role (catalyzing the conversion of methylcytosine to 5-hydroxymethylcytosine) and its established role in several aging/regenerative phenotypes [54, 55].

Several of the GWAS overlapping genes in the European ancestry analysis (e.g., TRIM59, KPNA4) are also implicated by our eQTL colocalization analysis. While it is important to note that these analyses do not test for causality, we identified regulatory relationships between SNPs associated with IEAA and expression of nearby genes. Moreover, four of these loci were not merely mQTLs for CpG sites that were used to construct the clocks. The IEAA-associated locus on chromosome 3 (lead SNP rs1047210) colocalized with an eQTL for TRIM59. DNA methylation levels at TRIM59 have been robustly associated with chronological age and its expression has been noted in multiple cancers [56-60]. The same locus was associated with IEAA in a previous GWAS (lead SNP rs11706810) [8]. In addition to TRIM59, the eQTL colocalization also implicated this locus in the expression of the pseudogene KRT8P12, along with SMC4 and KPNA4. There is evidence that SMC4 inhibits cellular senescence, an established hallmark of aging [61]. SMC4 is a core subunit of condensin complexes, which contribute to senescence processes, possibly through reorganization of genomic structure and transcriptional regulation [62]. KPNA4 is a member of the importin family of nuclear transport receptors, which work through the nuclear pore complex to selectively transport proteins from the cytoplasm to the nucleus. Like cellular senescence, dysfunction of nuclear transport has been proposed as a marker of aging [63]. Moreover, importin levels have been associated with impaired myocardial angiogenesis in aging [64]. A chromosome 1 locus associated with IEAA (lead SNP rs7550821) colocalized with an eQTL for CD46, encoding a regulator of T cell function and the complement system—a key component of the innate immune system where it promotes inflammation [65, 66]. There was also strong evidence for colocalization between an IEAA-associated locus on chromosome 15 (lead SNP rs12903325), and an eQTL for the lipid transporter gene ATP8B4, which contains variants that have been reported in relation to centenarian status in Italians and Alzheimer’s disease [67, 68]. An eQTL for CXXC4, encoding Idax, an inhibitor of Wnt signalling, colocalized with an IEAA-associated locus on chromosome 4 (lead SNP rs144317085). CXXC4 levels, along with KPNA4 and SMC4, have been associated with cancer progression [69-71].

Using the findings from the European ancestry meta-analyses, we observed shared genetic contributions between PhenoAge and GrimAge acceleration with education, cognitive ability, adiposity, and smoking. There were also significant epigenetic biomarker-specific genetic correlations with numerous other health and modifiable lifestyle factors. These are the first reported genetic correlations between GrimAge and PhenoAge acceleration and health and utilize much better powered GWAS results for IEAA and Hannum age acceleration. The best epigenetic biomarker of mortality risk, GrimAge acceleration, exhibited strong genetic correlations with parental longevity and lung cancer. As GrimAge uses a DNA methylation-based estimator of smoking pack-years in its definition, it is possible that the significant genetic correlation with lung cancer is mediated by smoking. Whereas many studies have observed phenotypic correlations between epigenetic clocks and health/aging outcomes, other than previous small-scale efforts [8, 9], none have done so through shared genetics. Furthermore, this is the first large-scale epigenetic clock study to attempt to separate correlation from causation. By comprehensively examining the genetic architectures of each clock, we begin to uncover the shared and unique biological signals that each clock is tracking.

Four different Mendelian randomization methods provide directional evidence of a causal influence of adiposity-related traits on GrimAge acceleration, while smoking cessation was inversely related to GrimAge acceleration. Several MR analyses indicate that increased educational attainment is associated with lower GrimAge acceleration. However, there was no causal evidence for associations with other lifestyle traits, such as alcohol consumption. There was limited evidence from MR to implicate any of the DNAm predictors as playing an important causative role in longevity or disease risk.

This is the first study to present polygenic risk scores for six epigenetic biomarkers of aging. A phenome-wide scan of the six polygenic risk scores (PRS) yielded highly significant associations between the PRS of GrimAge acceleration and adiposity-related traits, education, and parental longevity. Both the PRS analysis and the MR analyses suffer from two limitations: (i) the geographical structure in the UK Biobank cohort might confound these analyses [72] and (ii) low heritability estimates for some of the phenotypes (e.g., longevity or epigenetic biomarkers such as PAI1). The lack of replication of previous PAI1 GWAS findings, despite an equivalent sample size, questions the validity of the DNAm PAI1 proxy for GWAS analyses [25]. Limited statistical power due to low heritability or low sample size may help explain the disconnect between the genetic correlation analysis (which revealed a plethora of significant genetic correlations for GrimAge and PhenoAge) and the MR analysis (which led to a dearth of significant findings). In general, careful interpretation of the MR findings is required. Causal MR analyses that modelled epigenetic biomarkers as exposures and disease states as outcomes suffered from weak genetic instruments (e.g., for GrimAge acceleration and granulocyte proportions, where the variance explained was < 1%) or inadequate power in two-sample analysis. Future multivariate MR analyses will be required to test whether the protective causal effect of education on GrimAge is mediated by smoking, obesity, or other factors. These studies could also explore potential pleiotropy between the clocks, health outcomes, and white blood cell proportions. Furthermore, while the assumption of non-pleiotropic associations can be examined in the MR framework, this is not the case for the genetic correlation and polygenic risk score analyses.

Since carrying out these analyses, more accurate DNAm-based age predictors have been developed [73]. Compared to the original Hannum and Horvath clocks, the Zhang clock is less sensitive to cellular heterogeneity. Future studies should consider this clock for GWAS analysis. Furthermore, whereas we selected the original Hannum clock over a slightly modified version (labelled extrinsic epigenetic age acceleration), the correlation between the two age acceleration measures was > 0.95 in the two subsets (nSet1 = 5087; nSet2 = 4450) of the Generation Scotland study, the largest single cohort study in the meta-analysis.

Conclusions

Overall, this study highlights the shared genetic architecture between epigenetic aging, lifestyle factors (smoking, obesity), and parental longevity, which shows that DNAm-based biomarkers are valuable endophenotypes of biological aging.

Methods

Study cohort information

The meta-analysis sample comprised 36 datasets from 30 cohorts encompassing 40,905 participants (controls/healthy volunteers). Of these, 28 included individuals of European ancestry comprising 34,710 participants, and seven were of individuals of African American ancestry comprising 6195 participants. The total meta-analysis sample age range was 27.2–79.1 years (mean 54.0 years overall; 54.8 years across European ancestry cohorts; 50.4 years across African American cohorts) and comprised 0–100% females (mean 58.3%; 57.3% across European ancestry cohorts; 62.6% across African American cohorts). A Hispanic ancestry subset of the MESA cohort comprising 287 participants was used to compare effect sizes from trans-ethnic meta-analysis outputs and a European ancestry cohort comprising 1402 individuals (Young Finns Study) was used for polygenic prediction. Cohort-level descriptive data are presented in Additional file 2: Table S1 and described in Additional file 1.

Data preparation

Age-adjusted DNA methylation-based estimates of Hannum age, Intrinsic Horvath age, PhenoAge, GrimAge, plasminogen activator inhibitor-1 levels, and unadjusted granulocyte proportion were calculated using the Horvath epigenetic age calculator software (https://dnamage.genetics.ucla.edu/ or standalone scripts provided by Steve Horvath and Ake Lu). The following outputs were assessed: intrinsic epigenetic age acceleration—“IEAA”, Hannum age acceleration—“AgeAccelerationResidualHannum”, PhenoAge acceleration—“AgeAccelPheno”, GrimAge acceleration—“AgeAccelGrim”, estimate levels of Plasminogen Activation Inhibitor 1, adjusted for age—“DNAmPAIadjAge”, and estimated proportion of granulocytes—“Gran”. For each cohort, an outlier threshold for methylation values of +/− 5 standard deviations was applied and outlier samples were excluded from the analysis.

GWAS and meta-analysis

Quality control and imputation were done by each study separately (Additional file 1). Genotypes were imputed to either Haplotype Reference Consortium (HRC) or 1000 genomes phase 3 panels [74, 75]. In each cohort, association testing was conducted using imputed dosages using an additive model. Linear models were adjusted for sex and genetic principal components. GWAS summary statistics were obtained for between 1,097,816 and 15,221,271 genetic markers. This was the case for all cohorts with the exception of GOLDN (whole-genome sequence data) and the Sister Study (imputed data not available at the time of analysis). For each cohort, summary statistics were processed and harmonized using the R package EasyQC [76]. Multi-allelic variants were filtered to contain the variant with the highest minor allele count. At the individual cohort level, variants that were monomorphic, with a minor allele count ≤ 25, genotyped in < 30 individuals, or with an imputation quality score < 0.6 were removed. Allele codes and marker names were harmonized, duplicate variants were removed, and allele frequency checks were performed against the appropriate population reference data. Meta-analyses were performed with METAL using an inverse variance fixed effects scheme [77]. Meta-analyses were performed on European ancestry and African American studies separately (n = 34,710 and 6195, respectively). Variants were omitted from the meta-analysis if they were absent from > 50% of the total meta-analysis sample size. Cohort-specific genomic inflation factors ranged from 0.86 to 1.07. Genome-wide significance was defined as P < 5 × 10−8. To summarize the associations in terms of index SNP with the strongest association and other SNPs in linkage disequilibrium, we used conditional and joint association analysis of GWAS summary data, including the HLA region, in the GCTA-COJO software [22, 23]. A stepwise selection model was used with default settings for SNP LD (R2 < 0.9), analysis window size (10 Mb), and genome-wide significance (P < 5 × 10-8) using HRC imputed genotype data from Generation Scotland, and 1000G imputed genotype data from ARIC as the reference panels for the European ancestry and African American analyses, respectively. Heterogeneity I2 statistics were obtained from the meta-analyses and plotted against both −log10 p values and effect sizes to determine if SNPs with heterogeneous effects across cohorts were more statistically significant or had larger effect sizes. Systematic between-study heterogeneity was also investigated [21]. Meta-analyses were re-run after excluding cohorts identified as outliers and effect sizes were visually compared with the full meta-analysis output. Forest plots were prepared for all significant loci.

Trans-ethnic meta-analysis

A trans-ethnic meta-analysis of all European ancestry and African American cohorts was conducted using default settings in the Meta-Regression of Multi-Ethnic Genetic Association (MR-MEGA) tool [27]. We considered summary output for the first principal component of the meta-regression.

MR-MEGA summary statistics were uploaded to the Functional Annotation of Meta-Analysis Summary Statistics (FUMA) (http://fuma.ctglab.nl) software for annotation of the top loci using default settings, selecting 1000 Genomes phase 3 (all populations) as the reference population [79]. Independent lead SNPs had P < 5 × 10−8 and were independent of each other at r2 < 0.6; lead SNPs within this subset were required to have r2 < 0.1. A locus was defined by considering lead SNPs in a 250-kb range and all SNPs in LD (r2 ≥ 0.6) with at least one independent SNP.

Functional annotation of meta-analysis summary statistics

The European ancestry and African American meta-analysis summary statistics were uploaded to FUMA (http://fuma.ctglab.nl) for further annotation and functional analysis [78]. Genes were annotated from SNP-level data using the “SNP2GENE” tool, permitting gene set and tissue expression analyses using MAGMA [28]. Positional mapping was performed based on ANNOVAR annotations, applying a maximum distance of 10 kb between SNPs and genes. A Bonferroni-corrected significance threshold (adjusting for 18,606 tested genes) of P < 2.69 × 10−6 was set for the gene-based GWAS. Genes annotated to significant GWAS loci were further investigated using the “GENE2FUNC” tool in FUMA for enrichment of GWAS catalog gene sets, omitting the MHC region [79]. Bonferroni-corrected P value thresholds were applied.

Functional enrichment

To test if the GWAS meta-analysis findings were associated with regulatory and functional features of interest, enrichment analyses were conducted using GARFIELD [43]. SNPs were first pruned (r2 > 0.1) then annotated to categories (e.g., chromatin states, histone modifications, DNaseI hypersensitive sites, and transcription factor binding sites). Statistical enrichment was then carried out for SNPs at two P value thresholds (P < 1 × 10−5 and P < 1 × 10−8) while accounting for MAF, distance to the nearest TSS, and number of LD proxies.

Colocalization analysis

We hypothesized that some of the lead loci from the meta-analyses will have shared variants (1) across the Epigenetic biomarkers, (2) with DNAm sites in blood, and (3) with gene expression levels in blood. We used GoDMC summary statistics on 190,102 DNAm sites [29] to examine the overlap between loci and epigenetic clock DNAm sites, BMI-associated DNAm sites, and smoking-associated DNAm sites. We used eQTL Gen summary statistics on 19,942 transcripts [80] that were available in the MR-Base database [81]. We used the Rpackage gwasglue (https://mrcieu.github.io/gwasglue/) to extract SNPs that were +/− 1 Mb of the lead SNP and to harmonize the datasets. For each epigenetic biomarker–molecular trait pair (or pair of epigenetic biomarkers), we then performed colocalization analysis using the coloc.abf function in the R/coloc package [26], using default parameters. We only kept colocalized pairs with more than 50 shared SNPs and a posterior probability above 0.8 (PP > 0.80). We removed the HLA region from the eQTL colocalization analysis.

Disease and phenotype ontology enrichment

The potential role of Mendelian disease genes and associated pathways in influencing the epigenetic biomarkers was investigated with MendelVar [44], independently for each marker phenotype. We did not limit our analysis to any particular phenotype class among the Mendelian disease genes but looked for enrichment of any disease processes found to be strongly linked to genes in the GWAS loci. MendelVar analysis was run using intervals based on ± 0.5 Mbp window around the lead SNPs using the 1000 Genomes EUR population as LD reference [74]. Inside MendelVar, INRICH was run in “target” enrichment mode, with the target gene set filter set at minimum 5 (-i option) and maximum of 20,000 (-j option), and minimum observed threshold of 2 (-z option) [82]. The nominal P values were corrected for multiple testing with two rounds of permutation in INRICH.

Heritability and genetic correlation analysis

LD score regression, using LD scores and weights estimated from European ancestry populations (downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/), was used to assess genetic correlations between the six epigenetic biomarkers. Genetic correlations were further assessed between the six epigenetic biomarkers and publicly available GWAS summary statistics using the LDHub web interface (http://ldsc.broadinstitute.org/ldhub/) [45]. Meta-analysis results for each epigenetic biomarker were uploaded to the LDHub website, selecting all available traits for genetic correlation analysis. SNP heritability was estimated using univariate LD score regression. As the majority of large-scale GWA studies have been based on European ancestry populations, heritability and genetic correlation analyses were limited to this group to maximize statistical power. Filtering was performed to exclude traits where the LD Hub output came with the following warning messages: “Caution: using these data may yield less robust results due to minor departure of the LD structure” and “Caution: using this data may yield results outside bounds due to relative low Z score of the SNP heritability of the trait.” This left a total of 693 unique traits from 708 to 711 studies per epigenetic biomarker. A Bonferroni-corrected significance threshold of P < 0.05/693 = 7.21 × 10−5 was applied.

Polygenic risk scores

To determine the proportion of variance in the six epigenetic biomarkers that can be explained by common genetic variants, we carried out a polygenic risk score analysis using results from the European ancestry meta-analyses. Weights for the additive genetic scores were created by re-running the meta-analyses excluding one cohort (test cohort) at a time. Six weighted PGR scores (one for each epigenetic biomarker) were generated using default settings of the PRSice software (clump-kb = 250, clump-p = 1; clump-r2 = 0.25) [83]. P value thresholds were set at < 5 × 10−8, < 0.01, < 0.05, < 0.1, < 0.5, and 1. Linear regression models were built to calculate the incremental R2 between the null model (epigenetic biomarker~ sex) and full model (epigenetic biomarker ~ sex + polygenic risk score) in the test cohort. The procedure was iterated after excluding different test cohorts one by one (Lothian Birth Cohort 1921, Lothian Birth Cohort 1936, Framingham Heart Study, Born in Bradford, ARIES, and SABRE, respectively) from the meta-analysis. Finally, these steps were repeated, using the full meta-analysis summary statistic output to generate polygenic risk scores in a completely independent cohort (Young Finns Study, n = 1320).

For Born in Bradford, ARIES, and SABRE, best-guess genotype files with a MAF cut-off of 1% and info score > 0.8 were generated and the polygenic risk score analyses were corrected for 20 genetic PCs.

A phenome-wide association study of 581 heritable traits from the UK Biobank study was then carried out for polygenic risk scores based on independent SNPs with P < 5 × 10−8 or P < 1 from each of the six GWAS meta-analyses (http://mrcieu.mrsoftware.org/PRS_atlas/) [52]. The P value thresholds were based on the leave-one-out cohort PRS analyses described above (GrimAge acceleration and Granulocyte proportions: P < 1; IEAA, Hannum Age acceleration, PhenoAge acceleration, IEAA and PAI1 levels: P < 5 × 10−8). An FDR-corrected P value (PFDR < 0.05) was applied separately to each set of PheWAS results.

Mendelian randomization

To investigate if the epigenetic biomarkers were (i) causally influenced by lifestyle factors and (ii) had a causal effect on aging and disease outcomes, Mendelian randomization (MR) was performed in MR-Base [81]. The epigenetic measures were considered as both exposures (i.e., causally influencing the outcome) and outcomes (i.e., the epigenetic measure being causally influenced by a trait of interest). The analyses were further split into four sections: primary exposures/outcomes (common lifestyle risk factors and aging/disease outcomes from the largest available GWAS in MR-Base); secondary exposures/outcomes (traits identified as relevant via moderate genetic correlations from the LD regression analyses); 34 cell count exposures [24]; and 38 biomarker exposures [84]. SNPs instrumenting each exposure were clumped using a European LD reference panel and an r2 < 0.001. Harmonization of the SNP effects with the exposure and outcomes were performed so that palindromic SNPs were aligned when minor allele frequencies (MAFs) were < 0.3 or were otherwise excluded.

Inverse variance weighted (IVW) MR was carried out as the main analysis, with pleiotropy-robust sensitivity analyses featuring MR-Egger [85], weighted median [86], and weighted mode MR [87]. Significant associations were defined by a Bonferroni-corrected P value < 0.05. Where there was evidence for a causal effect based of the IVW model, we also assessed the potential for horizontal pleiotropy by means of heterogeneity assessment (Cochran’s Q-statistic) of individual SNP effects in both IVW and MR-Egger analyses and the Egger intercept test for directional pleiotropy [85].

Additional file 1. Individual cohort descriptions and acknowledgements.

Additional file 2. Supplementary Tables -Tables S1-S31.

Additional file 3. Supplementary Figures - Figures S1-S31.

Additional file 4. Assessment of genomic inflation and heterogeneity.

Additional file 5. Colocalization plots.

Additional file 6. Review history.