Epigenetic mutation load is weakly correlated with epigenetic age acceleration

Study population demographics

Our study includes 6,388 individuals from 4 studies: the Framingham Heart Study (FHS) Offspring Cohort, the Women’s Health Initiative (WHI), the Jackson Heart Study (JHS), and the Parkinson’s Environment and Genes (wave 1) known as the PEG1 study.

The main characteristics of the study populations are shown in Table 1. Briefly, FHS provided data for 2,326 individuals, with nearly half of them male (n = 1077; 46%) and all are white. Of the 2,091 female participants from the WHI, 989 (47%) are non-Hispanic white, 431 (21%) Hispanic, and 671 (32%) African American. JHS investigated 1,734 African American individuals with a majority of female participants (n = 1086; 63%). The 237 PEG1 control study participants were mostly non-Hispanic white (n = 207; 87%), and half were male (n = 126; 53%). The age ranges varied with the JHS having the largest range (22-93; mean = 56.2), and WHI the smallest (50-80; mean = 65.4). Mean ages of all populations ranged between 56.2 and 67.4. Additional details on cohorts and participant characteristics can be found in the Methods.

Epigenetic mutation load is the number of SEMs

All DNA methylation data was extracted from blood samples with the Illumina Infinium platform (450K array for PEG1, FHS, and WHI studies; EPIC array for WHI). Following a published and validated approach [24, 26, 31], a SEM is observed for a given person at a specific CpG site if an individual’s methylation level is more than three times the interquartile range (IQR) lower than the 25^th percentile (Q1 – 3 × IQR), or more than three times the IQR higher than the 75^th percentile (Q3 + 3 × IQR). The 25^th and 75^th percentile, and correspondingly the IQR, for each CpG locus was estimated across all samples. Furthermore, we defined the epigenetic mutation load (EML) of each study participant according to the total number of SEMs.

EML was highly variable across people (Supplementary Table 1), with a mean value ranging from 1647 to 3401 depending on the total number of CpGs measured on different arrays (FHS: 2433; WHI: 1647; JHS: 3401; PEG1: 2137). Since EMLs were not normally distributed, natural log-transformed EML values were used in all analyses.

EML was not associated with microarray slides (ANOVA p = 0.135) or position on the array (ANOVA p = 0.458). Also, EML was not correlated with the average intensity of bisulfite conversion controls (Pearson r = -0.085, p = 0.194). Thus, we concluded that the EML was independent of batches or other technical aspects.

Correlations among DNAm aging clocks

We calculated all DNAm aging estimators including the Horvath clock, the Hannum clock, the GrimAge clock, the PhenoAge clock, the SkinBlood clock, as well as an epigenetic estimate of telomere length (DNAmTL) using the online DNA Methylation Age Calculator (https://dnamage.genetics.ucla.edu/).

As expected, chronological age was strongly positively correlated with all epigenetic age estimators (Pearson r ranging from 0.79 to 0.93, Supplementary Figure 1), and these aging clocks were also strongly correlated with each other (Pearson r ranging from 0.73 to 0.90, Supplementary Figure 1). Meanwhile, the epigenetic estimate of telomere length, DNAmTL, was negatively correlated with chronological age and the epigenetic age estimates (Pearson r ranging from -0.63 to -0.72, Supplementary Figure 1).

For each clock, we calculated DNA methylation-based age acceleration based on the residuals of the regression of DNA methylation age on each participants’ chronological age. Thus, due to this approach, none of the epigenetic measures of age accelerations are correlated with chronological age (Pearson r = 0) as can be seen from Figure 1 and Supplementary Figure 2. AgeAccelHorvath is highly correlated (Pearson r = 0.93) with IEAA because both are based on the Horvath pan tissue clock. AgeAccelHannum was moderately associated with both AgeAccelHorvath and IEAA (Pearson r = 0.48 and 0.4 respectively). AgeAccelGrim showed only weak correlations with the other epigenetic measures of age acceleration which reflects the fact that GrimAge clock is a mortality risk predictor as opposed to an age estimator. Overall, the moderate pairwise correlations between the DNAm based biomarkers reflect different properties: some are highly confounded by blood cell composition and capture immunosenescence (Hannum, GrimAge, DNAmTL) while others are not (Horvath pan tissue, IEAA) [13, 32].

Figure 1. Heatmap of pairwise correlations of chronological age and epigenetic age accelerations. The heat map color-codes the pairwise Pearson correlations of chronological age and epigenetic age accelerations in the Framingham Heart Study (N=2326). Age represents the chronological age. AgeAccelHorvath, IEAA, AgeAccelHannum, and AgeAccelGrim represent measures of epigenetic age acceleration derived from the Horvath pan tissue clock, the intrinsic clock, the Hannum clock, and the GrimAge clock, respectively. The shades of color (blue, white, and red) visualize correlation values from -1 to 1. Each square reports a Pearson correlation coefficient.

Association between EML and DNAm aging clocks

We estimated the association between EML, chronological age, cell composition, and DNAm age acceleration using biweight midcorrelation (bicor) for each dataset separately and calculated pooled statistics using Stouffer’s method. Bicor is a median-based measurement of correlation that is robust to outliers [33]. We adjusted for potential confounders including age, sex, race/ethnicity, and cell compositions (naïve CD8 cells, CD8+CD28-CD45RA- T cells, Plasma Blasts, CD4 T cells, and Granulocytes) by regressing out the effects of these factors and retaining the residuals only for analysis. Results for AgeAccelHorvath, IEAA, AgeAccelHannum, and AgeAccelGrim are shown in Table 2 and Figure 2, while other clocks can be found in Supplementary Table 2. These analyses show that EML per study participant was positively correlated with chronological age (meta r = 0.171, meta P-value = 1.64E-42). Furthermore, EML was negatively correlated with CD4+ T cells (meta r = -0.121, meta P-value = 4.24E-22), plasmablasts (meta r = -0.085, meta P-value = 1.14E-11), and granulocytes (meta r = -0.064, meta P-value = 3.70E-07), but positively with exhausted CD8+ (defined as CD8+CD28-CD45RA-) T cells. These results are consistent with known age-related changes in blood cell composition [34, 35].

Figure 2. Correlations between EML and epigenetic age accelerations. Scatter plots of DNAm age acceleration estimators (x-axis; AgeAccelHorvath, IEAA, AgeAccelHannum, and AgeAccelGrim in each column, respectively) versus natural log-transformed EMLs (y-axis). Data from FHS, WHI, JHS, and PEG1 are plotted in four rows respectively. Each panel reports a biweight midcorrelation coefficient and correlation test p-value.

EML was also positively correlated with AgeAccelHorvath, IEAA, AgeAccelHannum, and AgeAccelGrim, with AgeAccelHannum exhibiting the strongest correlation (meta r = 0.179; meta P-value = 2.43E-46).

We further distinguished between epigenetic age acceleration and deceleration to determine correlations with EML. The correlation between EML and age acceleration was largely the same as what we presented originally. Interestingly, the correlation between EML and age deceleration was much smaller in size and less statistically significant (see Supplementary Table 3).

Sensitivity analyses

We evaluated associations between EML, chronological age, cell compositions, and age accelerations in males and females separately (Supplementary Table 4). For both sexes, EML remained positively correlated with chronological age, exhausted CD8+ T cells, and age acceleration suggesting that EML and age acceleration are independent of sex.

Several sensitivity analyses were conducted to ensure the reliability and reproducibility of the observed associations. To address a possibly non-linear relationship between epigenetic aging and chronological age, we additionally adjusted for a square term in age, (age², Supplementary Table 5). Also, to assess the potential for additional confounding, we adjusted for body mass index (Supplementary Table 6). This led to qualitatively similar results.

In order to explore whether the criteria used to define SEM will change results, we conducted another sensitivity analysis using two new SEM measures: 1) loose SEM: defined as a specific CpG site with its methylation level exceeding two times the interquartile range (IQR) of the first quartile (Q1 – 2 × IQR) or the third quartile (Q3 + 2 × IQR) across all subjects; and 2) stringent SEM: defined as a specific CpG site with its methylation level exceeding four times the interquartile range (IQR) of the first quartile (Q1 – 4 × IQR) or the third quartile (Q3 + 4 × IQR) across all subjects. We then calculated the total number of SEMs according to the loose and stringent definition for each person (loose or stringent EML, respectively). The biweight midcorrelations between loose or stringent EMLs and measures of epigenetic age accelerations were very similar to the original results (Supplementary Table 7).

We also explored the effect of different normalization methods for the methylation data (Illumina background correction, functional normalization, Noob, and quantile normalization). We found that the association between EML and age acceleration was not influenced by the normalization method (Supplementary Table 8).

Functional annotations

To test whether individual SEMs were randomly distributed across the genome or were more likely to be found in certain genomic regions or biological pathways, we conducted enrichment analyses to assess whether SEMs were enriched in clock-CpGs (Horvath clock, PhenoAge clock, Hannum clock), genomic regions (transcription start sites (TSS1500, TSS200), untranslated regions (5’UTR, 3’UTR), 1^st Exon, and gene body), or regulatory features (i.e., enhancers, DNase hypersensitive sites, open chromatin regions, transcription factor binding site, promoters). For each participant, we first annotated the probes and each mutation based on the location related to genes (i.e., TSS1500, TSS200, 5’UTR, 1^st EXON, gene body, 3’UTR), or regulatory features using the manifest provided by Illumina. Then we conducted hypergeometric tests for each region and each subject separately using a nominal significance threshold of 0.05. Last, for each region, we summarized the number of individuals for which the test was significant (Supplementary Tables 9–11).

The results show that for each clock-CpG set or region, only a small proportion of participants have their SEMs enriched which illustrating both the stochastic nature and inter-personal variation of SEMs.

Pathway enrichment analysis

Next, we examined whether, among study participants who exhibit faster age accelerations, SEMs are enriched in particular biological pathways. We first conducted KEGG pathway enrichment analyses for each study participant. Then for each KEGG pathway, we calculated the number of people with enriched SEMs for this particular pathway. Finally, we investigated the association between SEMs enrichment in a particular pathway and age accelerations using linear regression. Additional details can be found in the method section.

Supplementary Tables 12–16 showed the top 10 pathways that were commonly enriched in each study, and generally we found that SEMs enriched in these pathways were also statistically significantly associated with faster age accelerations (AgeAccelHorvath, IEAA, AgeAccelHannum, and AgeAccelGrim, respectively).

Briefly, in all four study populations, we identified similar enrichment patterns, and SEMs enriched in signaling pathways, axon guidance, glutamatergic synapse, morphine addiction, glucocorticoid pathway (Cushing syndrome), or circadian rhythm pathways were associated with faster AgeAccelHorvath, AgeAccelHannum, and IEAA. Whereas the associations between pathway enrichment and AgeAccelGrim or AgeAccelPheno were less strong and not necessarily statistically significant. We only observed SEMs enriched in neuroactive ligand-receptor interaction was associated with faster AceAccelGrim and AgeAccelPheno.

Region-specific EML

To address the functionality of SEMs on biological age acceleration, we calculated the number of SEMs co-located with clock CpGs for each study participant (i.e., clock-specific SEMs) and assessed whether there were any clock-specific EMLs corresponding to age acceleration (Supplementary Table 17), but we observed no statistically significant association with the three clocks tested (Horvath clock: 353 CpGs; Hannum clock: 71 CpGs; PhenoAge clocks: 513 CpGs).

Next, we divided CpGs into different genomics region/regulatory feature groups based on the annotations, and then calculated EMLs within each region for each study participant (i.e., genomic region-specific EML; regulatory region-specific EML) (Table 3 and Supplementary Table 18). EMLs in TSS1500, TSS200, and the 1stExon regions were related to faster age accelerations. Also, EML in DNase hypersensitive regions was positively correlated with faster age accelerations. In contrast, 3’UTR specific-EML was associated with younger chronologic age and slower age acceleration.

The direction of SEM

Based on the direction of the mutation, we separated SEMs into hypomethylated SEM (Q1 – 3 × IQR) and hypermethylated SEM (Q3 + 3 × IQR) and calculated hypomethylated and hypermethylated EMLs respectively within FHS. We assessed the correlations between the newly calculated directional EMLs and epigenetic age acceleration. The results remained largely the same (See Supplementary Table 19). Furthermore, consistent with previous studies [29], hypermethylated SEMs were mainly located in CpG islands, while hypomethylated SEMs were enriched in the open seas (see Supplementary Tables 20, 21).

Shannon entropy, EML, and DNAm age acceleration

As an alternate measure for a well-functioning EMS that maintains genomic stability, we calculated the Shannon entropy of the whole methylome based on the 450K or EPIC array. An increase in entropy means that the methylome becomes less predictable across the population of cells, i.e. when the methylation fractions (beta values) tend towards 50%.

We found the Shannon entropy to be positively correlated with chronologic age (meta r = 0.046, meta P-value = 2.16E-04), EML (meta r = 0.234, meta P-value = 7.71E-78), and all four measures of age accelerations (meta P-value: AgeAccelHorvath = 8.79E-09, IEAA = 6.56E-04, AgeAccelHannum = 1.80E-22, AgeAccelGrim = 1.67E-22) (Table 4 and Supplementary Table 22).

Study population

Our study is based on data from four studies: the Framingham Heart Study (FHS) Offspring Cohort, the Women’s Health Initiative (WHI), the Jackson Heart Study (JHS), and the Parkinson’s Environment and Genes (wave 1) known as the PEG1 study.

We used 2,326 individuals from the FHS Offspring cohort [48]. The FHS cohort is a large-scale longitudinal study started in 1948, initially investigating the common factors of characteristics that contribute to cardiovascular disease (CVD) (https://www.framinghamheartstudy.org/index.php). The study at first enrolled participants living in the town of Framingham, Massachusetts, who were free of overt symptoms of CVD, heart attack, or stroke at enrollment. In 1971, the study started FHS Offspring Cohort to enroll a second generation of the original participants’ adult children and their spouses (n= 5124) for conducting similar examinations. The FHS Offspring Cohort collected medical history and measurement data, immunoassays at exam 7, and blood DNA methylation profiling at exam 8. Participants from the FHS Offspring Cohort were eligible for our study if they attended both the seventh and eighth examination cycles and consented to have their molecular data used for the study. We used the 2,326 participants from the group of Health/Medical/Biomedical (IRB, MDS) consent and available for both Immunoassay array DNA methylation array data. The FHS data are available in dbGaP (accession number: phs000363.v16.p10 and phs000724.v2.p9).

The WHI is a national study that enrolled postmenopausal women aged 50-79 years into the clinical trials (CT) or observational study (OS) cohorts between 1993 and 1998 [49, 50]. We included 2,091 WHI participants with available phenotype and DNA methylation array data from “Broad Agency Award 23” (WHI BA23). WHI BA23 focuses on identifying miRNA and genomic biomarkers of coronary heart disease (CHD), integrating the biomarkers into diagnostic and prognostic predictors of CHD and other related phenotypes, and other objectives can be found in https://www.whi.org/researchers/data/WHIStudies/StudySites/BA23/Pages/home.aspx.

The JHS is a large, population-based observational study evaluating the etiology of cardiovascular, renal, and respiratory diseases among African Americans residing in the three counties (Hinds, Madison, and Rankin) that make up the Jackson, Mississippi metropolitan area [51]. The age at enrollment for the unrelated cohort was 35-84 years; the family cohort included related individuals >21 years old. Participants provided extensive medical and social history, had an array of physical and biochemical measurements and diagnostic procedures, and provided genomic DNA during a baseline examination (2000-2004) and two follow-up examinations (2005-2008 and 2009-2012). Annual follow-up interviews and cohort surveillance are ongoing. In our analysis, we used the visits at baseline from 1,734 individuals as part of project JHS ancillary study ASN0104, available with both phenotype and DNA methylation array data.

The PEG1 study was conducted during 2000-2007 to investigate the causes of Parkinson's disease (PD) in agricultural regions of the California central valley. We analyzed blood samples from 238 healthy controls enrolled from Kern, Tulare, or Fresno counties. Controls were required to be over the age of 35, having lived within one of the counties for at least 5 years prior to enrollment, and do not have a diagnosis of Parkinsonism. Demographic information, lifestyle factors, and medication use were collected in standardized interviews, including lifetime information of cigarette smoking and coffee/ tea consumption.

For PEG1, FHS, and WHI, peripheral blood samples were collected, and bisulfite conversion using the Zymo EZ DNA Methylation Kit (Zymo Research, Orange, CA, USA) as well as subsequent hybridization of the HumanMethylation450k Bead Chip (Illumina, San Diego, CA), and scanning (iScan, Illumina) were performed according to the manufacturer’s protocols by applying standard settings. For JHS, DNA methylation quantification was conducted using HumanMethylation EPIC Bead Chip (Illumina, San Diego, CA).

Preprocess

For PEG1 samples, raw signal intensities were retrieved using the function read.metharray.exp of the R package minfi from the Bioconductor open-source software (http://www.bioconductor.org/), followed by linear dye bias correction, noob background correction, and functional normalization using the same R package [52–55]; β-value was used for all the analyses. One sample was identified as low quality due to low median methylated and unmethylated signal intensities across the entire array and thus removed from the study population. Detection p-values were derived using the function detectionP as the probability of the total signal (methylation + unmethylated) being detected above the background signal level, as estimated from negative-control probes. All in all, 845 probes with a detection p-value above 0.05 in at least 5% of samples were removed. Also, 645 probes with a bead count <3 in at least 5% of samples; 11,334 probes on the X or Y chromosome; 7,306 probes containing a SNP at the CpG interrogation site and/or at the single nucleotide extension for 5% maf; and 27,332 cross-reactive probes were also removed. In total, 438,050 probes were included for downstream analyses.

For FHS and WHI samples, 11,334 probes on the X or Y chromosome; 7,306 probes containing a SNP at the CpG interrogation; and 27,332 cross-reactive probes were also removed. In total, 439,540 probes were included for downstream analyses.

For JHS samples, 19,532 probes on the X or Y chromosome; 53,435 probes containing a SNP at the CpG interrogation and cross-reactive were also removed. In total, 793,869 probes were included for downstream analyses.

SEM calculation

The calculation of SEM was consistent with a previously published and validated approach [24, 26, 31]. CpG with methylation levels three times the interquartile range above the third quartile or below the first quartile was identified as a SEM. Toward this end, we calculated the IQR for each of the 438,050 probes in each dataset (for PEG1, FHS, and WHI) or the 793,969 probes (for JHS). Then, SEMs were identified based on extreme methylation levels. Finally, we summed across the count of all SEMs per sample and defined the total number of SEMs of each study participant as epigenetic mutation load (EML). EML was not normally distributed; therefore, we used the natural log of the number of SEMs for all regression analyses.

In FHS, we separated SEMs into hypermethylated and hypomethylated SEMs based on the direction of the mutation. We also defined consistently hypermethylated or hypomethylated SEMs as a CpG mutated in the same direction in more than 10 participants.

In order to assess whether the criteria used for SEMs will change the results, we defined loose SEM and stringent SEM as described above. We then calculated the total numbers of loose and stringent SEMs for each person.

DNA methylation age

We included eight different DNAm aging biomarkers in this study. Utilizing our online DNA Methylation Age Calculator (https://dnamage.genetics.ucla.edu/), we calculated DNA methylation-based ages and the age accelerations based on the residuals of the regression of DNA methylation age on each participants’ chronological age for each clock.

Four types of DNA methylation-based biomarkers were included in the main analyses. Briefly, Horvath clock was calculated using a linear combination of 353 CpGs that have previously been shown to predict chronological age in multiple tissues [5]; and the intrinsic clock was derived from the Horvath clock by additionally regressing out cell compositions [30]; Hannum clock was calculated using a linear combination of 71 CpGs to predict chronological age in blood [11]; and GrimAge clock) was calculated from a linear combination of 7 DNAm plasma protein surrogates and a DNAm-based estimator of smoking pack-years designed to predict mortality [7].

Other DNAm aging biomarkers were included in the Supplementary analyses, including: the extrinsic clock [30], PhenoAge clock [6], SkinBlood clock [56], DNAm based estimator of telomere length [57], each of the 7 DNAm protein surrogates underlying the definition of the GrimAge clock [7], as well as DNAm based estimate of smoking pack-years.

Cell composition

White blood cell composition was imputed for each study participant using our online published DNA Methylation Age Calculator, https://dnamage.genetics.ucla.edu/. The following imputed blood cell counts were included in downstream analyses: CD4+ T, naïve CD8+ T, exhausted cytotoxic CD8+ T cells (defined as CD8 positive CD28 negative CD45R negative, CD8+CD28-CD45RA-), plasmablasts, and granulocytes. Naïve CD8+ T, exhausted cytotoxic CD8+ T cells, and plasmablasts were calculated based on the Horvath method [58]. The remaining cell types were imputed using the Houseman method [59].

Shannon entropy

The formula of Shannon entropy is:

$\begin{array}{l}Entropy=\frac{1}{N\times lo{g}_2(\frac{1}{2})}\\ \times {\displaystyle \sum _{i=1}^N[{\beta }_i}\times lo{g}_2({\beta }_i)+(1-{\beta }_i)\\ \times lo{g}_2(1-{\beta }_i)]\end{array}$

where β_i represents the methylation beta value for the i^th probe (CpG site) in the array, N represents the total number of probes included in the formula [11].

Statistical analysis

All analyses were conducted using R v3.6.1. We used Pearson correlations to assess the relations between different DNAm ages and age accelerations. We evaluated potential batch effects by assessing the difference of EMLs in microarray slides or position on the array with the ANOVA test. To eliminate the possibility that SEMs are driven by incomplete bisulfite conversion, Pearson correlations between EMLs and the average intensity of bisulfite conversion controls were also calculated. Average intensity of bisulfite conversion controls was derived using the ENmix R package [60].

To assess the association between EML and age/age accelerations/cell compositions, we applied biweight midcorrelation (bicor) implemented in the WGCNA R package. We adjusted for potential confounders including age, sex, race/ethnicity, and cell compositions (naïve CD8 cells, CD8+CD28-CD45RA- T cells, Plasma Blasts, CD4 T cells, and Granulocytes) by regressing out the effects of these factors and retaining the residuals of log(EML)only for analysis.

We also conducted stratified analyses to evaluate the associations between EML, chronological age, epigenetic estimates of cell composition, and epigenetic age acceleration in males and females separately for the PEG1, FHS, and JHS studies. To ensure the reliability and reproducibility of the associations, several sensitivity analyses were conducted. In addition to the potential confounders mentioned above, we adjusted for a quadratic term in age (age^2) to account for non-linear relationships and body mass index.

To evaluate the effect of data preprocessing steps, we repeated the analysis using four different normalization methods (Illumina background correction, noob, functional normalization, and quantile normalization implemented in the minfi package) in the PEG dataset.

We analyzed each dataset separately, therefore in order to obtain an overall p-value across all four studies, we conducted meta-analyses using Stouffer’s method for meta-analysis estimates of the correlation coefficient (meta r), and the corresponding two-sided p-values (meta p-value).

Functional annotation, region-specific SEMs, and pathway enrichment analysis

For each participant, we annotated the probes and SEMs based on the location in relation to genes (TSS1500, TSS200, 5’UTR, 1^st Exon, gene body, and 3’UTR), or regulatory features (enhancer region, DNase hypersensitive region, open chromatin region, transcription factor binding site, and promoter region) using the manifest provided by Illumina. To test for region specific enrichment, we conducted hypergeometric tests for each region and each subject separately. A p-value less than 0.05 was considered statistically significant. For biological pathway enrichment analysis, the enrichKEGG function implemented in the ClusterProfiler package [61] was used to assess whether study participant’s SEMs were enriched in particular KEGG pathways (P-value threshold = 0.05). For each genomic region or KEGG pathway, we then summarized how many study participants had their SEMs enriched. We also investigated the association between SEMs enrichment in each KEGG pathway and age accelerations using linear regression, adjusted for the total number of SEMs per study participant:

$\begin{array}{l}Forpathwayj:AgeAcce{l}_i\\ ={\beta }_0+{\beta }_{1}Enric{h}_{ij}\\ +{\beta }_2\log (TotalEM{L}_i)+{\epsilon }_i\end{array}$

where AgeAccel stands for age accelerations (AgeAccelHorvath, IEAA, AgeAccelHannum, AgeAccelGrim); Enrich stands for enrichment for pathway j (significant: 1, non-significant: 0); Total EML stands for log transformed total EML for participant i.