DNA methylation clocks for estimating biological age in Chinese cohorts

To establish an epigenetic clock that would enable investigations of underlying associations with multi-omics data, we collected whole blood DNA methylation data from a group of individuals aged between 20 and 87 years old (Quzhou cohort from southern China, 112 women, 138 men, Table S1), and subjected biological samples from this cohort to phenomics, transcriptomics, proteomics, metabolomics, and facial image analysis (Li et al., 2023a) (Fig. 1A). In addition, we used a larger cohort (CAS cohort from northern China, 639 women, 431 men) as a validation of the established DNA methylation clock (Du et al., 2019) (Fig. 1A). This comprehensive dataset allowed us to thoroughly profile age-related DNA methylation alterations and build epigenetic predictions into multi-modal aging clocks in the Chinese population.

In the initial principal component analysis (PCA), we identified a significant association between age and DNA methylation (Fig. 1B). We also observed a significant reduction in global DNA methylation level along with an elevation in entropy of methylation in the elderly of both genders (Fig. 1C and 1D). Taken together, these results suggest aging is associated with dramatic changes in DNA methylation. Then, we calculated the age-related differential methylation positions (DMPs) (change in beta per year of age > 0.002 and BH-adjusted P-value < 0.05, corrected with blood cell proportions, gender, and BMI values) to identify 331 and 1,330 DMPs that were upregulated and downregulated with age, respectively (Fig. 1E; Table S2). After annotating these age-related DMPs, we observed that DMPs upregulated with age were primarily located within CpG island regions, whereas DMPs downregulated with age were predominantly isolated CpGs situated outside any CpG island, falling into the open sea region (Fig. 1E). Additionally, the percentage of hypermethylated sites situated within 200 bp upstream of the transcription start site (TSS) reached 16% (Fig. 1E). These findings indicate that hypermethylated sites may be more likely to regulate gene expression during aging (Marttila et al., 2015; Pal and Taylor, 2016). Moreover, pathway enrichment analysis suggested that DMPs upregulated with age are involved in various biological processes, particularly pathways associated with development, as identified in previous research (Lu et al., 2023a), while those downregulated with age are enriched for cell morphogenesis (Fig. S1C).

Of note, we observed a comparatively higher concentration of DMPs in the TSS regions of chromosomes 3, 6, and 19 compared to other chromosomes (Fig. 1H and 1I). Within chromosome 3, the hypermethylated sites were primarily situated in the regions of the GRM2 gene (Fig. 1I), whose methylation level has been reported positive associated with age (Naue et al., 2017). In chromosome 6, the hypomethylated sites were predominantly located in the regions of the ZC3H12D gene (Fig. 1I), which encodes MCPIP4, a CCCH-zinc finger family protein known to regulate the pro-inflammatory activation of macrophages (Liang et al., 2008). The hypermethylated sites of chromosome 6 were mainly found in the gene segments of ELOVL2, which has been strongly associated with aging in various tissues (Garagnani et al., 2012a; Slieker et al., 2018). In chromosome 19, the methylation level of NACHT and WD repeat domain containing 1 (NWD1) was observed to be lower in the elderly (Fig. 1I), and its DNA methylation levels have previously been suggested as a marker for aging (Chen et al., 2022). Collectively, these findings highlight the genomic sites on chromosomes 3, 6, and 19 as hotspots for age-related DNA methylation changes.

In addition to linear change with age, we also tried to identify crests of DNA methylation alterations across age stages. In sliding window analysis, we found that the number of DMPs peaked at the age of 35 in both genders (Fig. 1F; Table S3). Furthermore, women had an additional peak of methylation changes at age 55 that was not observed in men (Fig. S2F). Although the number of DMPs at 55 is nearly twice the number of DMPs at 35 in women, the number of hypomethylated DMPs was high in both age stages (Fig. S2G), suggesting that loss of DNA methylation might occur at earlier stages of lifespan (Figs. 1G,S2F, S2G, S3F, and S3G).

To interpret the biological implication of the aging epigenome, we performed an integrative analysis of DNA methylation and both bulk and single-cell transcriptomic data from the same cohort. This strategy allowed us to identify age-related probe-gene pairs (Yao et al., 2015), including key transcription factors (TFs) potentially whose activities regulated by DNA methylation, and estimate how these genes change with age in a cell type-specific manner (Fig. 2A). In total, we obtained 207 age-upregulated DMPs with downregulated nearby genes (127 in total), and 53 age-downregulated DMPs with upregulated nearby genes (39 in total) (Fig. 2B; Table S4). In gene set variation analysis (GSVA) based on the scRNA-seq data, we demonstrate that these 166 genes were preferentially changed in T-cell subsets (e.g., Naïve/Memory CD8T and Naïve CD4T) (Fig. 2C; Table S5), suggesting that the DNA methylation state of T cell populations is sensitive to aging. By comparing these genes with genes documented in the Aging Altas (Aging Altas, 2021) we found that several known aging-related genes were potentially regulated by the altered DNA methylation (Fig. 2E). For example, CDKN2A, a canonical marker for cellular aging (Krishnamurthy et al., 2004), is upregulated during aging and hypomethylated in the elderly (Fig. 2E), consistent with the results recorded in the Lineage Landscape (Yan et al., 2023). Furthermore, the activated probe-gene pairs (decreased DNA methylation and increased expression) were enriched in pathways such as platelet activation, Ras protein signal transduction, while the repressed probe-gene pairs (increased DNA methylation and decreased expression) are associated with glycogen biosynthetic process and positive regulation of cold-induced thermogenesis, both of which were closely related to homeostasis (Ngo and Steyn, 2015; Wang et al., 2022) (Fig. 2D).

To predict which TFs might bind to genomic loci with aging DMPs, we examined the overlap between probe-gene pairs identified here and TFs whose recognition sequence motifs were documented in the HOCOMOCO database (Kulakovskiy et al., 2016) (see Methods). The results inferred that a larger number of TFs were enriched on upregulated DMPs than on downregulated DMPs (Fig. 2F; Table S6). TFs predicted to bind upregulated DMPs with the highest enrichment scores included NRF1, ZBT14, KAISO, MECP2, etc. (Fig. 2F). Methyl CpG binding protein 2 (MeCP2) acts as a transcriptional regulator by binding to methylated DNA on CpG islands to repress transcription (Meehan et al., 1992), supporting that our analysis can identify potential TFs contributing to the regulatory network of the aging methylome. Consistently, loss of nuclear respiratory factor 1 (NRF1), a key regulator in proteostasis maintenance (Cui et al., 2021; Lehrbach and Ruvkun, 2019), is thought to contribute to aging especially in naïve T cells, causing loss of chromatin accessibility at gene promoters (Moskowitz et al., 2017a, 2017b). We also noticed that several repressed genes in the identified probe-gene pairs were NRF1 target genes, and associated with proteostasis, including NPM1, NMD3, USP444 (Matsuura-Suzuki et al., 2022; Thakur et al., 2015; Zhang et al., 2012) (Figs. 2G andS1E). In all, these results indicate that NRF1 regulates the loss of proteostasis in a DNA methylation-dependent manner.

We also screened for TFs for which expression levels were also changed with aging, and in an inverse direction with DNA methylation level of their binding sites. We only identified five TFs that met this criterion: AHR, E2F6, INSM1, PATZ1, and PLAG1 (Figs. 2G andS1G). Interestingly, the aryl hydrocarbon receptor (AHR) was identified as a stress-induced DNA methylation reader (Habano et al., 2022) associated with the CD8T cells activation (Zhang et al., 2023c). Additionally, E2F6 and ZBT14 co-regulate transcriptional expression with AHR, and are also enriched in the hypermethylated DMPs (Oshchepkova et al., 2020; Puga et al., 2000). These results imply that AHR and its co-regulators play a key role in age-related methylation changes. Taken together, our findings establish a DNA methylation-dependent TF regulatory network that is associated with various aging-related pathways and provides insights into the epigenetics mechanism of immune cell aging.

We then built a DNA methylation clock (iCAS-DNAmAge) based on the Quzhou cohort and applied the model performance to healthy individuals of the CAS cohort to validate its robustness (see Methods) (Fig. 3A; Table S7). The iCAS-DNAmAge comprises 65 CpG sites, with 35 being upregulated with age and associated with protein ubiquitination, a pathway closely related to protein homeostasis that is critical in the aging process (Huang et al., 2023; Zhang et al., 2022), while the remaining 30 are downregulated with aging and linked to protein phosphorylation (Fig. 3D and 3E). The results showed that prediction accuracy is high in both the Quzhou and CAS cohorts (R_Quzhou = 0.97, MAE_Quzhou = 3.45, R_CAS = 0.77, MAE_CAS = 4.37) (Fig. 3B). To test the generalization ability of the model, we also compared our model to previous DNA methylation clocks (Hannum et al., 2013; Horvath, 2013; Levine et al., 2018). We were encouraged to find that the performance of the iCAS-DNAmAge model surpassed those of the three previous clocks in both Quzhou (R_{Hannum’s clock} = 0.97, MAE_{Hannum’s clock} = 11.20; R_{Horvath’s clock} = 0.96, MAE_{Horvath’s clock} = 10.38; R_PhenoAge = 0.95, MAE_PhenoAge = 6.26) and CAS cohort (R_{Hannum’s clock} = 0.79, MAE_{Hannum’s clock} = 11.69; R_{Horvath’s clock} = 0.76, MAE_{Horvath’s clock} = 9.92; R_PhenoAge = 0.75, MAE_PhenoAge = 7.61) (Fig. 3B and 3C; Tables S8 and S9). This result demonstrates that the iCAS-DNAmAge model is generalized in different Chinese cohorts and suggests that, although the DNAm clock can usually predict age accurately, ethnicity might also affect its age predictions.

To gain mechanism insights into the iCAS-DNAmAge clock, we compared the components of all four models and only found seven probes shared between the iCAS-DNAmAge and previous clocks (Fig. 3F). The methylation levels of probes located at ELOVL2, FHL2, KLF14, NHLRC1, and TRIM59 increase with age, while that at MKLN1 decreases with age (Fig. 3F). The CpG sites located at FHL2 appeared in all four clocks above, and their DNA methylation levels have been frequently reported to be associated with aging (Bacalini et al., 2017; Garagnani et al., 2012b; Jung et al., 2019), indicating their close relationship with the aging process. Interestingly, a recent study reported that the combination of methylation levels of ELOVL2, FHL2, KLF14, TRIM59, and C1orf132/MIR29B2C is sufficient for age prediction (Jung et al., 2019). Based on these findings, we picked the 7 CpG sites incorporated in iCAS-DNAmAge and the previous clock, developed a concise DNAm clock with a stepwise regression strategy (Fig. 3G). Ultimately, by using only five CpG sites, including cg16867657 (ELOVL2), cg22454769 (FHL2), cg08097417 (KLF14) (Table S7), we obtained an age prediction model (5CpG-DNAmAge) that achieved an accurate prediction in both the validation set of Quzhou cohort (R = 0.95, MAE = 3.90) and the CAS cohort (R = 0.82, MAE = 5.37) (Fig. 3H; Table S10).

Furthermore, the genes related to the CpG sites underlying iCAS-DNAmAge were predominantly enriched in export from the cell, RHO GTPase cycle, and protein phosphorylation pathways (Fig. 3E). Among them, RHO GTPase has been shown to impact the proliferation and differentiation of hemocytes and is strongly associated with the senescence of human lymphoblastoid cell lines (Florian et al., 2017; Mulloy et al., 2010). Additionally, protein phosphorylation has been found to affect blood cell morphology and cell adhesion (D’Alessandro et al., 2015), implying a close link between iCAS-DNAmAge and the aging process.

Previously, we generated a group of multi-modal aging clocks for female individuals from Quzhou cohort, constructed based on multi-omics and facial image datasets (Li et al., 2023a). Named after their features, these clocks encompass compositeAge, facialAge, transAge, hormoneAge, proteinAge, phenoAge, metabAge, immuneAge, and lipidAge. For instance, compositeAge is derived from the combination of phenomics, transcriptomics, proteomics, and metabolomics datasets, facialAge is generated based on analysis of facial features, and immuneAge is constructed utilizing immune-related features from the transcriptome, proteome, and phenotype. Next, we asked whether these clocks, which reflect aging from both distinct and composite perspectives, could be replaced by the DNA methylation data, which represents an easier method for evaluating biological age by simply estimating the DNA methylation of blood cells (Fig. 4A). Using the leave-one-out cross-validation strategy, we predicted multi-modal ages for female individuals in the validation set (see Methods). Encouragingly, the DNA methylation was proved to be a robust surrogate for these aging clocks. For compositeAge, facialAge, and transAge specifically, which have a high correlation with chronological age, we found that their corresponding DNA methylation clocks were correspondingly more accurate (Fig. 4B and 4C). For example, compositeAge has the best prediction power of chronological age over other clocks based on a single data type, and the prediction results of compositeAge-DNAmAge, also achieve the highest accuracy (R = 0.92, MAE = 3.35) (Fig. 4B and 4C). Similarly, facialAge-DNAmAge (R = 0.88, MAE = 4.34), hormoneAge-DNAmAge (R = 0.88, MAE = 4.65), phenoAge-DNAmAge (R = 0.88, MAE = 4.90), transaAge-DNAmAge (R = 0.85, MAE = 4.33), proteinAge-DNAmAge (R = 0.85, MAE = 4.71), and metabAge-DNAmAge (R = 0.82, MAE = 5.13) models exhibit strong correlations with their respective age indices. Though associated with lower accuracy, models built from DNA methylation data can also predict the other clocks, including immuneAge-DNAmAge (R = 0.72, MAE = 6.49), lipidAge-DNAmAge (R = 0.53, MAE = 7.98) (Fig. 4B and 4C; Table S9).

Subsequently, we applied the DNAm multi-modal clocks to the CAS cohort. We found that the prediction results of compositeAge-DNAmAge and proteinAge-DNAmAge were well correlated with the chronological age (Fig. 4B and 4C; Table S9). The same high correlations were also identified between predicted multi-modal ages and chronological age in the Quzhou cohort (Fig. 4B and 4C). These findings indicate that DNA methylation-based clocks could present a streamlined and consistent method for proxying the complex data provided by multi-omic clocks.

The methylation sites located on FHL2, TRIM59, ELVOL2, and MKLN1 appeared in multiple methylation clocks, further suggesting a strong correlation between these sites and chronological age (Fig. S1F). Pathway enrichment analysis of clock-specific sites showed that genes related to these sites were enriched in pathways including regulation of Wnt signaling pathway, hematopoiesis, etc., suggesting the DNAm multi-modal clocks may have the ability to reflect different levels of aging (Fig. S1H).

Having established the DNAm multi-modal clocks, we postulated that these DNA methylation-based clocks can shed light into what factors contribute to the heterogenous aging process. When we estimated the age pace for different DNAm multi-modal clocks, we did not observe any significant correlation between the age pace and chronological age (Figs. 5B, S5A; Tables S11, S12, S13 and S14), suggesting that individual variations in biological aging were not influenced by age (Elliott et al., 2021). We then identified individuals in the two tails (20%) of the age pace distribution as age-decelerator and age-accelerator. By comparing the multi-omics data, especially the clinical performances of these two groups, we were able to identify key health indicators related to the pace of different aging clocks (Fig. 5A). The results indicated that the compositeAge-DNAmAge exhibited the highest number of multi-omics features associated with its age pace, followed by transAge-DNAmAge and facialAge-DNAmAge (Fig. S6A), suggesting their strong correlation with health status.

Elevated chronological inflammation is an important hallmark of aging (López-Otín et al., 2023; Wu et al., 2024). Here, based on several aging clocks, we found that higher values in age-accelerators were associated with inflammation and immune activation. For example, inflammation-related proteins were upregulated in age-accelerators defined by compositeAge-DNAmAge (CSF1, EXT1, ITGB2), facialAge-DNAmAge (BST1, DPP4), immuneAge-DNAmAge (FCN1, ERAP1), hormoneAge-DNAmAge (IL6R) as well as iCAS-DNAmAge (IL6R, HRNR) (Fig. S6B; Table S12), suggesting the important role of immunosenescence in aging process (Liu et al., 2023; Zhao et al., 2023). Of note, higher cytomegalovirus (CMV) antibody levels, were reported to be associated with various aging-related diseases (Merani et al., 2017; Parry et al., 2016) and accelerated aging in previous studies (Kananen et al., 2015; Poloni et al., 2022; Stowe et al., 2012), was correlated with faster aging based on almost all aging clocks, including compositeAge-DNAmAge, transAge-DNAmAge, metabAge-DNAmAge, and hormoneAge-DNAmAge (Figs. 5B andS6C), implying its instrumental role in affecting the rate of aging. We also found that several proteins down-regulated in the facialAge-DNAmAge accelerated group, including CSTA, KRT2, and TGM3, were associated with epidermal cell differentiation and skin development (Fig. S6B). These results denote how DNAm-based clocks can indicate divergent aging-related molecular changes that link to critical pathological process, such as activated inflammation and impaired immune cell function.

At the phenotypic level, blood-related parameters: hemoglobin (Hb), mean corpuscular hemoglobin (MCH), and mean corpuscular hemoglobin concentration (MCHC) displayed elevation in compositeAge-DNAmAge accelerators (Figs. 5B andS6C). Elevated levels of these parameters have been reported to be associated with chronic hypoxia, like chronic obstructive pulmonary disease and heart disease (Martínez-Quintana and Rodríguez-González, 2013; Woodson et al., 1970). Additionally, transAge-DNAmAge accelerators exhibited elevated MCH and MCHC, along with elongated QRS interval (Figs. 5C and S6B). In facialAge-DNAmAge accelerated group, extracellular fluid/total body fluid (ECF/TBF) and extracellular water/total body water (ECW/TBW) ratios were increased, which were positively correlated with mortality (Kim et al., 2017; Pérez-Morales et al., 2021) (Figs. 5B andS6C). Our analysis revealed that lifestyle choices contributing to an age-accelerating effect were characterized by reduced sleep duration in the context of immuneAge-DNAmAge, lipidAge-DNAmAge, and metabAge-DNAmAge (Fig. 5C). Additionally, for compositeAge-DNAmAge, a diet comprising daily intake of pickled and fried foods was identified (Fig. 5C). Notably, within the extensive CAS cohort, we observed that individuals with hypertension exhibited a significantly higher pace of compositeAge-DNAmAge compared to healthy individuals (Fig. 5D). Taken together, these data demonstrate that DNAm-based clocks can serve as a more convenient alternative to estimate biological age and predict health indicators, and might even be applied to monitor one’s disease status (Fig. 5E).

The Quzhou cohort was conducted in Quzhou, Zhejiang Province of China. All volunteers signed a written informed consent. Volunteer information was obtained from a previous study (Li et al., 2023a).

A 15 mL whole blood sample was collected from each volunteer after fasting overnight. Two EDTA-covered anticoagulation tubes (for plasma and blood cells) were used for collection and transported immediately to the laboratory and stored at 4°C pending further processing. All blood samples were stored at −80°C for long-term preservation.

DNA was obtained with the TIANamp Genomic DNA Kit (Tiangen, DP304-03). The quality of the DNA was determined by checking the OD₂₆₀/OD₂₈₀ ratio in the spectrophotometer and the integrity in agarose gels, which should be between 1.8 and 2.0, and the length of the DNA should be greater than 10 kb.

Genomic DNA (gDNA) samples were bisulfite converted with the EZ DNA Methylation™ Kit (Zymo Research, USA) and then hybridized to the Infinium Methylation EPIC BeadChip array (Illumina, USA) and the Infinium Methylation EPIC v2.0 BeadChip array (Illumina, USA) according to the manufacturer’s instructions. Cases and controls were randomly assigned to each array. The BeadChip array was run in a single base extension reaction, stained, and imaged on an Illumina iScan.

DNA methylation data for each sample was first processed with ChAMP (version 2.21.1) (Tian et al., 2017), to perform out-of-band Infinium I intensities background correction, RELIC dye bias correction, inter-array normalization, RCP probe type bias correction, low-quality loci filtering, and imputation with the k-nearest neighbor (KNN) method. The beta value matrix was normalized with the BMIQ method (Marabita et al., 2013), and the batch effect was removed with champ.runCombat function (Johnson et al., 2007).

SNP-related probes, multi-hit probes, and probes on the X and Y chromosomes were then deleted referring to previous studies (Nordlund et al., 2013; Zhou et al., 2017), and only the overlapping probes of the two chips were preserved for the pan-gender analysis. All remaining CpG sites were annotated based on the EPICanno.ilm10b4.hg reference (Hansen, 2017).

The probes on the sex chromosomes were retained for methylation analysis of males and females separately (with the X chromosome for females and both the X and Y chromosomes for males). This process involved using the EPICanno.ilm10b4.hg reference annotated for female CpG sites and the manifest file from the Illumina official website for male CpG sites.

CpG sites were annotated based on their position on the genome and CpG island. In this context, CpG island refers to a short stretch of palindromic DNA with a “CpG” sequence. Additionally, CpG shore represents a region within 2 kb around the CpG island, CpG shelf refers to a region within 2 kb around the CpG shore, and Opensea denotes a region outside the aforementioned regions. Furthermore, TSS 200 denotes the region 200 bp upstream of the transcription start site, while TSS 1500 represents the region 1,500 bp upstream of the transcription start site. Moreover, 5ʹUTR and 3ʹUTR represent the 5ʹ and 3ʹ untranslated sequences, 1stExon denotes the first exon, ExonBnd represents the exonic region, Body signifies the gene body region, and IGR stands for the intergenic region (Hansen, 2017).

The beta value for each methylated site was calculated using the formula [M/(M+U)], where M represents the methylation intensity of the site, and U represents the unmethylated intensity of the site, in the range of 0 to 1.

ENmix R package (version 1.36.03) (Kassebaum et al., 2016) was used to estimate the proportions of CD8⁺ T cells, CD4⁺ T cells, NK cells, monocytes, B cells, and neutrophils based on FlowSorted.Blood.EPIC dataset (Salas et al., 2018) with the houseman method (Houseman et al., 2012).

The lmFit function of the limma R package (version 3.54.2) (Ritchie et al., 2015) was applied to fit linear models for the beta value of each methylation site with age as a continuous variable and the proportion of blood cells and BMI as covariates. The P-value was corrected using the Benjamin-Hochberg method. The CpG sites with BH-adjusted P-value < 0.01 and |ΔBeta/year| > 0.002 were identified as age-related DMPs.The genomicDensity function from the circlize package (version 0.4.15) (Gu et al., 2014) was used to identify the distribution of DMPs within 1M bp bins.

The DEswan package (version 0.0.0.9001) (Lehallier et al., 2019) was used to conduct sliding window analysis for detecting peaks in age-related methylation changes. Statistics analysis were conducted with a sliding window size of 5 years in parcel of 5 years. Those with BH-adjusted P-value < 0.05 were considered as significantly changed methylation sites.

Genes in probe-gene pairs were selected to perform GSVA with the GSVA R package (version 1.5.0) (Hänzelmann et al., 2013). GSVA scores for different cell types between aged and young individuals were calculated based on the scRNA-seq datasets generated from the Quzhou cohort (Li et al., 2023a). Differences in GSVA scores between young and aged groups were considered significant when P-value < 0.05 (t-test).

The Metascape webtool (Zhou et al., 2019) was used to conduct pathway enrichment analysis for genes and proteins.

Statistical analysis of the comparisons in Figs. 2C, 5C, 5D, S4A, S4B and S6B was performed using the two-tailed t-test with the stats R package (version 4.2.2). Statistical analysis of the comparisons in Fig. 5E was performed using the Wilcoxon rank-sum test with the ggpubr R package (version 0.4.0).

DNA methylation raw data of Quzhou cohort can be accessed in the OMIX database under the accession number OMIX005455. DNA methylation raw data of CAS cohort will be accessed with reasonable requests. The multi-omics datasets were obtained from a previous study (Li et al., 2023a).