Integrative epigenetics and transcriptomics identify aging genes in human blood

Domain: biology
Published: 2026-01-19
Editorial outcome: major_revision

Abstract
Recent epigenome-wide studies have identified a large number of genomic regions that consistently exhibit changes in their methylation status with aging across diverse populations, but the functional consequences of these changes are largely unknown. On the other hand, transcriptomic changes are more easily interpreted than epigenetic alterations, but previously identified age-related gene expression changes have shown limited replicability across populations. Here, we develop an approach that leverages high-resolution multi-omic data for an integrative analysis of epigenetic and transcriptomic age-related changes and identify genomic regions associated with both epigenetic and transcriptomic age-dependent changes in blood. Our results show that these multi-omic aging genes in blood are enriched for adaptive immune functions, replicate more robustly across diverse populations and are more strongly associated with aging-related outcomes compared to the genes identified using epigenetic or transcriptomic data alone. These multi-omic aging genes may serve as targets for epigenetic editing to facilitate cellular rejuvenation.

Introduction
Aging is associated with a wide range of molecular and cellular changes. Understanding the impact of these changes may provide insights into the biology of aging and help identify potential therapeutic strategies for aging and age-related diseases. DNA methylation (DNAm) alterations rank among the most well-studied and prominent age-related molecular changes; importantly, these alterations generally replicate well across many different cohorts and species1. DNAm is part of the epigenome and is involved in regulating gene expression, but the downstream consequences of most age-related DNAm changes remain incompletely understood, rendering it challenging to interpret their functional impact. Conversely, transcriptomic changes are more functionally informative than DNAm due to their indication of the current level of gene expression in tissues and cells. Age-related changes to the transcriptome have also been described and have facilitated the development of ‘transcriptomic clocks’2; however, these clocks generally feature lower replicability across cohorts and datasets than DNAm-based clocks3, higher noise, and susceptibility to batch effects4. These features make comparison of existing transcriptomic clocks across studies difficult.

The integration of multiple omics-based measures could assist in overcoming specific limitations of an individual readout. Previous studies have leveraged multi-omics analyses to derive DNAm-based predictors (or proxies) of features in other omics or clinical data, combining the reliability of DNAm profiling with the more functional insights derived from other molecular changes. For instance, DNAm-based EpiScores have been shown to predict plasma proteomic features5, and recently, a DNAm-based multi-omic aging clock was described6. The combination of epigenetic and transcriptomic data has been suggested to improve biological age predictions2, yet few studies to date have investigated direct links between age-related epigenetic and transcriptomic changes1,7–11. These studies have reported conflicting findings, with some reporting that methylome changes are indicative of alterations in gene expression7–9 and others noting limited association with expression of affected genes1,10,11.

A comprehensive, integrative analysis of epigenetic and transcriptomic changes across the human lifespan has not yet been conducted. Such a study could shed further light on the functional relevance of age-related DNAm changes and identify genes that could be targeted to combat aging-related declines in physiological health. Here, we perform a large-scale analysis of DNAm and gene expression data from blood samples, utilizing data from several cohorts (n = 4174 and 3461 total samples for DNAm and RNA-seq, respectively; Table 1). Our integrative multi-omic approach enabled the identification of functionally relevant DNAm changes in blood associated with gene expression alterations, hereafter referred to as multi-omic aging genes. We validate our findings in an external high-quality blood DNAm dataset, generated using the latest methylation array technology, from the Mass General Brigham Biobank (n = 500) that comprises a broad range of ages and has additional data on many aging-associated outcomes. Our findings provide a first in-depth assessment of the functional consequences of age-related DNAm changes in blood, and highlight the potential for multi-omics approaches to uncover functionally relevant genes and genetic loci. In turn, these features could be used to develop meaningful predictors for relevant aging-related outcomes and may be targeted by therapeutics to mitigate aging-related molecular changes and decline.

Results
Age-associated gene expression is less reproducible than DNAm
We initially aimed to identify age-associated genes from six large transcriptomic datasets, including data from the Multi-Ethnic Study of Atherosclerosis (MESA), which features two sampling timepoints (referred to hereafter as MESA1 and MESA2); the Parkinson’s Progression Markers Initiative (PPMI); the Gates Grand Challenge (GC6); the 500 Functional Genomics Project (500FG); and the JenAge Ageing Factor Database. These datasets cover wide age ranges and are balanced between male and female participants (Fig. 1a). For each cohort, we assessed the correlation of gene expression levels with age (see “Methods”). Among the top 10 transcripts associated with age across all datasets, we observed that the correlation with age was quite variable between cohorts (Fig. 1b). For instance, we found loss of expression of CD248 as a function of age, consistent with previous reports12,13, but the magnitude of the correlation ranged from −0.3 to −0.5 (Fig. 1c). For further validation of this observation, we examined the correlation between age and expression levels of aging transcripts previously identified by Peters et al. across these cohorts. In agreement with our previous result, we found generally low levels of correlation between expression levels of these genes and age across cohorts (Fig. 1d). Thus, age-related transcriptomic data appear to replicate poorly across cohorts.

We next asked whether the same was true of DNAm data. We analyzed DNAm from six large data sources: MESA1, MESA2, and PPMI, as above; a cohort of patients with Rheumatoid Arthritis and healthy controls (RA); the Grady Trauma Project (Grady); and the Genetic Epidemiology Network of Arteriopathy (GENOA) study (Fig. 1e). For the top 10 CpG sites that were found to be associated with age, we found considerably higher and more consistent correlations with age across the datasets (Fig. 1f). For instance, cg16867657, a well-known aging-associated CpG site within the ELOVL2 promoter14, emerged as the top site consistently positively correlated with age (Fig. 1g). Across aging CpGs reported by Varshavsky et al. in 202315, we again found correlations that were larger in magnitude than those observed for aging transcripts, and more consistent across cohorts. Thus, age-associated DNAm data are considerably more replicable across cohorts than transcriptomic data. Statistical comparison of correlation strengths confirmed that age-associated DNAm data show significantly higher and more consistent correlations with age across datasets compared to aging transcripts (Supplementary Fig. 1).

Integrating omics identifies validated age-associated genes
To identify age-associated genes that are consistent across blood samples of several cohorts, we integrated DNAm and transcriptomic data from studies that produced both data modalities (MESA and PPMI) to identify genes whose age-associated DNAm changes have functional consequences at the level of gene expression. On average, across MESA timepoints and across cohorts, we observed surprisingly little correlation between gene expression levels and age for genes whose CpG sites had the greatest gain of methylation over aging (Fig. 2a), consistent with some previous reports1,10. Conversely, we found that genes with the greatest decrease in expression levels with aging consistently featured significantly higher DNAm correlation with age (Fig. 2b). This observation was statistically significant across the MESA1, MESA2, and PPMI cohorts (Supplementary Fig. 2). These observations were consistent at the level of individual genes: There were minimal differences in the expression of ELOVL2, KLF14, and FHL2, the genes associated with the top-ranked aging-associated CpGs, between young and old individuals (Fig. 2c, left, see Methods for details). In contrast, the top down-regulated aging transcripts, CD248, LRRN3, and NELL2 each featured increased DNAm in older compared to younger individuals (Fig. 2c, right). CD248 in particular showed strong positive correlation with age at the level of DNAm (Fig. 2d, left) and strong inverse correlation with gene expression (Fig. 2d, right); this was true at both MESA time points and in the PPMI dataset. Interestingly, we found our multi-omic aging genes to be particularly enriched for T cell-specific genes (e.g., CD27, CD28, CD248, TCF7).

Gold-standard reference DNAm dataset for aging blood
To validate our identified multi-omic aging genes, we aimed to perform an independent external validation using an external cohort. To demonstrate that our integrative analyses and the identified epigenetic-transcriptomic aging genes are replicable, we performed an independent external validation using an independent cohort we specifically designed for omic analysis of aging biomarkers across the lifespan. We generated DNAm profiles for 500 individuals of diverse ages from the Mass General Brigham (MGB) Biobank using the Illumina Infinium Methylation EPIC v2.0 array, which covers over 935,000 CpG sites enriched for regulatory regions and has been shown to exhibit high reproducibility16. Our subjects were recruited from a major metropolitan academic medical center, were roughly balanced between male and female, and were generally representative of the racial/ethnic distribution of the local area (Fig. 3a).

As expected, DNAm levels at aging CpGs15 exhibited consistent age-dependent changes, with slightly stronger correlations for the MGB cohort compared to other cohorts (Fig. 3b). The age-CpG associations in the MGB 500 cohort were robust, showing high concordance between Pearson and Spearman correlation methods, indicating that results are not driven by outliers (Supplementary Fig. 3). Sites with the strongest DNAm correlation with age, particularly those that gain methylation, were enriched in lowly-expressed and repressed genes (Fig. 3c). Such genes were enriched for strong PRC2 binding (Fig. 3d, e), consistent with our previous analyses of genome-wide PRC2-binding and genomic regions with age-associated CpGs17. Further enrichment analyses of these PRC2-associated multi-omic aging genes highlighted pathways related to immune function, cell differentiation, and developmental processes (Supplementary Fig. 4). We also confirmed the gain of methylation at CpG sites within promoters of genes whose expression is inversely associated with age in this dataset (Fig. 3f). Finally, we again compared the concordance of the top age-associated genes identified by Peters et al. with gene expression data, with the addition of DNAm data from the MGB cohort. We were able to validate 106 genes that feature both an age-dependent gain of DNAm and loss of expression levels; we term such genes multi-omic aging genes (see Methods, Supplementary Data 1). For these 106 multi-omic aging genes, there was no strong overall association between the magnitude of age-correlation for gene expression and the age-correlation for their associated CpG methylation (Supplementary Fig. 5). The addition of DNAm data considerably increased the percentage of age-associated genes we were able to validate in both MESA timepoints and in PPMI, as well as in the other three cohorts with RNA-seq data (Fig. 3g). KEGG pathway enrichment analysis confirmed significant enrichment for T-cell specific functions among these multi-omic aging genes, including T cell receptor signaling and differentiation pathways (Supplementary Fig. 6). In contrast, the broader set of 1497 age-associated genes identified by Peters et al. based on transcriptomics alone did not show particular enrichment for T cell-specific pathways (Supplementary Fig. 7). To further explore the regulatory landscape, transcription factor enrichment analysis using ChEA identified key upstream regulators, including PRC2 components (e.g., EZH2, SUZ12) and factors involved in immune/haematopoietic cell development (e.g., RUNX), coordinating the expression of these multi-omic aging genes (Supplementary Fig. 8). These analyses demonstrate that integration of multiple omics modalities has the ability to refine the identification and validation of true aging-associated genes, where both DNA methylation and gene expression are jointly altered, pointing towards functional effects. Moreover, the high-quality dataset we generated exhibits strong age-related DNAm changes, covers a broad age range, and features a multitude of aging-associated outcome data such as mortality. Thus, this dataset is ideal for benchmarking existing and future biomarkers of aging.

Multi-omic aging genes predict aging outcomes.
To evaluate the functional relevance of the identified multi-omic aging genes, we analyzed their associations with aging-related outcomes, particularly mortality. We utilized large datasets containing both DNA methylation and mortality data, including a subset of participants from the MGB Biobank (MGB-4K) and the Generation Scotland (GS) cohort. CpG sites associated with our multi-omic aging genes were strongly associated with all-cause mortality risk in both cohorts (Fig. 4a), with hazard ratios for the top CpGs ranging from approximately 1.3–1.75 per standard deviation (Fig. 4b). Notably, these CpGs selected via our multi-omic approach showed stronger and more consistent associations with mortality compared to CpGs selected based solely on age correlation without transcriptomic filtering (Supplementary Fig. 9). Importantly, survival analysis based on the top mortality-risk-associated CpG site showed clear stratification of subjects (Fig. 4c). Kaplan-Meier plots for additional top CpG sites in both the MGB-4K and GS cohorts, with and without cell composition adjustment, are provided in Supplementary Fig. 10 and Supplementary Fig. 11, respectively. In addition, disease gene network over-representation analysis of our multi-omic aging genes revealed enrichment for numerous aging-associated diseases, particularly those related to the aging immune system, such as various lymphomas (Fig. 5f). Thus, multi-omic aging genes are strongly associated with clinically relevant aging-related outcomes such as disease and mortality risk.

Age-related methylation in aging genes is independent of cell composition changes
Blood is a heterogeneous tissue composed of multiple cell types, the proportions of which change with age. Given these age-dependent alterations in blood cell composition18, age-related changes observed in the transcriptome and epigenome may potentially be related to changes in cell proportions over time19. Noting a clear enrichment for T cell differentiation genes in our multi-omic aging genes (Fig. 5a, b), we sought to investigate whether aging patterns of our multi-omic aging genes were truly intrinsic to individual blood cell types or simply a reflection of age-related cell composition changes. To do so, we generated methylation-based cell composition estimates20,21 for our MGB cohort samples (Fig. 3a) using predictors of proportions for the most abundant blood cell types and ordinal abundance values for less common, yet age-dependent, cell types. We observed modest age-related changes in cell composition for the most common cell types (Fig. 5c). However, for naive T cells, we found a dramatic decline associated with age that is consistent with a well-documented feature of human immune aging (Fig. 5d)22. Finally, we characterized the relationship between age and methylation level for CpGs within our identified multi-omic aging genes after adjusting for cell composition. Significant associations of methylation levels at these CpGs with age remained after the adjustment (Fig. 5e), indicating that age-related changes in these CpGs are independent of changes in cell composition. The results were consistent when adjusting for cell type compositions using 12 cell-type deconvolution (Supplementary Fig. 12). Thus, methylation levels at CpG sites within our multi-omic aging genes are robustly associated with aging.

Methods
All research was conducted in accordance with relevant ethical regulations. The study protocol was approved by the Mass General Brigham Institutional Review Board (protocol 2021P003059) and the Tayside Committee on Medical Research Ethics, on behalf of the National Health Service, for Generation Scotland.

Study datasets
The current study leverages several previously-described DNA methylation and RNA sequencing datasets, shown in Table 1. Their underlying cohorts or population characteristics have been described in detail in publications by the original respective authors. In addition, we generate a DNA methylation dataset encompassing 500 individuals of various ages derived from the Mass General Brigham Biobank. This study was approved by the Mass General Brigham IRB (protocol 2021P003059). Details on methylation data and RNA sequencing data processing are described below. For comparison of molecular data between young and old groups, participants 45 years old or younger are labeled as young, and those 75 years or older are labeled as old.

The Generation Scotland (GS) cohort comprises 18,859 participants (58.8% women, 41.2% men) aged 18–99 years, recruited between 2006 and 2011. Whole blood samples were collected, DNA extracted, and methylation measured on the Illumina Infinium Methylation EPIC v1.0 array. Mortality follow-up was ascertained through the National Records of Scotland up to October 2023 (median follow-up 13.9 years). Data were processed using standard quality control pipelines, including background correction, probe filtering, and normalization25.

The Mass General Brigham 4K (MGB4K) cohort included 4243 participants (65.4% women, 34.6% men). Whole blood samples were collected, DNA was extracted, and methylation profiles were generated using the Illumina Infinium Methylation EPIC v1.0 array. Mortality follow-up was ascertained through linkage with electronic health records until July 28, 2022 (10.8% deceased). Data preprocessing included background correction, probe filtering, and normalization using the R SeSAMe package6.

MGB500 samples were selected from healthy control donors of diverse ages from the MGB Biobank. The donors were recruited from a major metropolitan academic medical center, were roughly balanced between male and female, and were generally representative of the racial/ethnic distribution of the local area (Fig. 3a).

MESA cohort is a population-based sample of an initial 6814 men and women aged 45–84 without known cardiovascular disease. Thirty-eight percent of the recruited participants were White, 28 percent African American, 22 percent Hispanic, and 12 percent Asian, predominantly of Chinese descent. Participants were recruited from six field centers across the United States: Wake Forest University, Columbia University, Johns Hopkins University, University of Minnesota, Northwestern University and University of California - Los Angeles. (mesa-nhlbi.org).

Public DNAm and RNA-seq data were selected among commonly used datasets used in epigenetic or transcriptomic studies of aging from the generally healthy populations, with the exception of PPMI, which was included to increase our multi-omic data. Diseased participants from the GC6 cohort were excluded from the analyses.

Participant characteristics and consent
This study included human participants from multiple cohorts. No participants were recruited for this study. Blood samples from MGB Biobank were used to generate DNA Methylation data. The Mass General Brigham (MGB500) cohort comprised 500 individuals (461 after quality control) between 20 and 90 years of age, with a roughly balanced distribution of males and females. The MGB4K cohort consisted of 4243 individuals aged 20–90 years, of whom 65.4% were female, and 34.6% were male. Generation Scotland included 18,859 individuals aged 18–99 years, of whom 58.8% were female, and 41.2% were male. The Multi-Ethnic Study of Atherosclerosis (MESA) recruited 6814 participants aged 45–84 years, including White, African American, Hispanic, and Asian individuals with a balanced sex distribution. Additional cohorts, including PPMI, GC6, 500FG, JenAge, RA, Grady, and GENOA, are described in Table 1.

Participants provided written informed consent at the time of recruitment, as required by each biobank or study protocol. Ethical approval was obtained from the Mass General Brigham Institutional Review Board (protocol 2021P003059) and from the Tayside Committee on Medical Research Ethics, on behalf of the National Health Service, for Generation Scotland. Participants were not financially compensated for contributing biospecimens to the MGB Biobank or the Generation Scotland resource, in accordance with institutional policies for biobank participation.

Sex of participants was determined by self-report at enrollment, and gender identity data were not collected. Sex information was available and analyzed in aggregate when reporting cohort demographics. The study was not specifically designed to test sex- or gender-related differences in methylation or gene expression, since the primary aim was to identify aging-associated changes across populations, and therefore analyses were not stratified by sex. Future work will be required to investigate sex-specific effects in more detail. This study complies with SAGER guidelines: sex was reported for all cohorts, and justification has been provided for not conducting sex-specific analyses.

Human tissue samples
Human peripheral blood samples for DNA methylation profiling were obtained from the Mass General Brigham (MGB) Biobank (MGB500 and MGB4K cohorts). All participants provided informed consent for the research use of their biospecimens and data. Samples were processed to generate DNA methylation profiles using the Illumina Human Methylation EPIC arrays. Residual DNA samples have been conserved within the MGB Biobank and may be available for further analysis. Requests for access to these samples can be made to the Mass General Brigham Biobank (biobank@partners.org) and will be considered in accordance with institutional and IRB policies. Samples from Generation Scotland (GS) were obtained under approval of the Tayside Committee on Medical Research Ethics, and have been retained in the Generation Scotland resource (contact: access@genscot.ac.uk). No unique samples were destroyed for this study beyond routine laboratory consumption during DNA extraction and array processing.

Methylation data generation and processing
Previously processed described datasets were obtained from relevant sources (e.g., GEO, dbGaP, or others), and no further processing or filtering was conducted. Methylation data from the MGB biobank were generated using the Illumina Human Methylation EPIC version 2 array (cat# 20087709), encompassing over 950.000 CpG sites. Raw .IDAT files were preprocessed using standard parameters in the R SeSAMe package26, version 1.22.1, using the recommended processing for EPICv2 data, with the QCDPB parameter (qualityMask, inferInfiniumIChannel, dyeBiasNL, pOOBAH, noob). Samples and CpGs with more than 10 percent missing values were excluded from the generated beta value matrix. PRC2 binding regions were identified using EZH2 ChIP-seq peaks from ENCODE (ENCFF109KCQ), with CpGs within these regions marked as PRC2-associated sites, as in Moqri et al. 2024. Quality control analysis examined the distribution of the plate number as well as Type I and Type II probe loadings on the leading principal component across all cohorts to verify the absence of systematic technical bias (Supplementary Fig. 13).

RNA seq data and processing
We obtained previous RNA-seq datasets as raw gene expression counts. Count data were normalized using log transformation on (raw counts+1), and the median adjusted to account for batch effects when merging datasets. The majority of analyses was conducted within individual datasets rather than merging multiple datasets.

Integrating DNA Methylation and RNA-seq data to identify multi-omic aging genes
To integrate DNA Methylation and RNA-seq data, we employed a 4-step process:
1) We initially identified aging transcripts (ATs) in the cell type of interest, as defined by a statistically significantly Pearson correlation with chronological age (p-value < 0.05), starting from the 1497 transcripts described as the ageing transcriptome by Peters et al. (2015);
2) We outlined aging regulatory regions (ARGs) for each AT identified in step 1. ARGs were defined as the 1500 bases flanking the transcription start site (TSS) in Gencode Hg38;
3) Aging CpGs (ACs) within each ARG were identified by marking regions with n ≥ 3 consecutive cytosine-guanine (CpG) sites whose DNA methylation levels were significantly associated with chronological age at R2 > 0.2;
4) Finally, candidate aging genes (candidate AGs) were defined by filtering ACs whose methylation levels exhibited a discordant directionality with their associated ATs, e.g., a gain in DNAm with an associated downregulation in expression with age.

This process resulted in 106 multi-omic aging genes (Supplementary Data 1).

Validated aging genes (validated AGs) were defined as candidate AGs if the direction of age-dependent changes was confirmed in independent datasets. Of note, the number of validated genes in each cohort varies.

Survival analysis
Linkage to records for the MGB-4k was described previously in Chen et al. 20236. GS survival data was ascertained via linkage to the National Records of Scotland. The dates of death range from the time of the baseline assessment (that took place, depending on participant, between 2006 and 2011) and October 2023, with the latter date representing the date of censoring for alive participants. The median follow-up for GS data is 13.9 years (alive participants: 14.0; dead participants: 9.4).

For both datasets, we applied a multivariate Cox Proportional Hazard regression model to test the association between all-cause mortality and each aging CpG site, adjusting for age and sex. Next, we ranked the aging CpG sites by p-values from smallest to largest and selected the top 10 sites, as demonstrated in the volcano plot. For these top aging CpGs, we drew a forest plot to show the point estimates and 95% confidence intervals of the adjusted hazard ratios for all-cause mortality. In addition, we plotted the adjusted survival curves for these top 10 aging CpGs.

Cell composition analysis
To estimate the cell composition of each sample in the MGB cohort, we applied software provided by the Clock Foundation that deconvolutes cell composition using bulk methylation profiles. This software estimates the proportion of cell types through constrained quadratic programming for common blood cell types (CD8+ T cells, CD4+ T cells, NK cells, B cells, monocytes, and granulocytes)21 and penalized regression for rarer cell types (plasmablasts, CD8+ CD28- CD45RA- T cells, naive CD8+ T cells, and naive CD4+ T cells)20.

For each CpG site associated with a multi-omic aging gene, we corrected for cell composition through the following approach: first, we applied multivariate linear regression to predict the methylation beta value of that CpG site using, as inputs, age and the cell proportion of each cell type, excluding granulocytes and plasmablasts due to high collinearity (defined as a variance inflation factor over 5). Next, we determined whether the regression coefficient p-value for age was still significant (Benjamini-Hochberg-corrected p-value < 0.05) in predicting CpG methylation.

Functional enrichment analysis. Functional enrichment of genes associated with aging across modalities was performed using Fisher’s exact test and gene terms from GO Biological Process, using all available genes in the ontology with the clusterProfiler package in R. Redundant functions with high semantic similarity were filtered out with the simplify function. Adjustment for multiple comparisons was performed with the default Benjamini-Hochberg approach, and gene terms with adjusted p-value < 0.05 were considered statistically significant. Enriched functions were visualized with a dot plot. Gene-concept network of multi-omic blood biomarkers of aging and representative enriched function terms was constructed with the cnetplot function from the enrichplot package in R. Over-representation analysis for the disease gene network (Fig. 5f) was conducted using R version 4.3.3 and the enrichDGN function from the DOSE package. This function leverages DisGeNET (Janet et al., 2015) to construct disease-associated gene networks.

Statistics & reproducibility
No statistical method was used to predetermine sample size. Sample sizes were determined by the availability of participants in each cohort and are consistent with previous large-scale epigenetic and transcriptomic studies. No data were excluded from the analyses unless they failed predefined quality control criteria (e.g., > 10% missing CpGs or poor sequencing quality), as described above. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

All statistical analyses were conducted in R (version 4.3.3) using standard packages as detailed above. Associations between DNA methylation or gene expression and age were assessed using Pearson correlation and validated with Spearman correlation for robustness. Multiple testing correction was performed using the Benjamini–Hochberg method, with adjusted p-values < 0.05 considered significant. Survival analyses were conducted using multivariate Cox proportional hazards regression, adjusting for age and sex.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In \"Integrative epigenetics and transcriptomics identify aging genes in human blood,\" Moqri et al. present an integrative analysis of age-associated DNA methylation and RNA expression signatures. The authors can pinpoint roughly 100 genes robustly associated with aging by analyzing DNA methylation and RNA expression data from various cohorts. Importantly, enrichment analyses suggest that the identified genes may provide a better starting point for studying the molecular mechanisms of aging in the human blood. The study achieves this goal by generating an independent methylation data set to validate findings from previous, partially publicly available studies. Naturally, the generation of an independent RNA expression data set of the same cohort would have been an added benefit - although not strictly required, in my opinion, to substantiate the conclusions and claims of this study.\nConceptually, the manuscript introduces the term \"multi-omic aging genes.\" While it may not be entirely new to seek genes exhibiting both transcriptional and regulatory changes during processes in time, the manuscript makes an essential case for analyzing multiple aging-associated features simultaneously. It demonstrates that this approach helps to identify biologically meaningful \"aging genes.\"\nThe manuscript is well-written, and the presented data substantiates all conclusions and interpretations. Nevertheless, I have some suggestions that the authors may find helpful in revising their manuscript.\n- The authors may like to discuss the advantages and limitations of DNAm arrays. For instance, arrays can measure only a fraction of all genomic CpGs, thus substantially limiting the detection of age-associated DNAm alterations. On the other hand, arrays may exhibit higher reproducibility, e.g., compared to RRBS or WGBS.\n- The current manuscript implicitly goes beyond integrating DNAm and RNA expression already. Specifically, the authors show the enrichment of age-associated CpGs in regions affected by PRC2 binding. An increasing number of studies show that DNAm changes occurring during aging or in age-associated diseases are often associated with PRC2-associated silencing. Also, because PRC2 plays a role in bivalently marked genes that frequently play a role in development, I'd be curious whether the subset of multi-omic aging genes with PRC2 association enriches any more specific functions.\n- The integration of DNAm and RNA expression data is straightforward but somewhat conservative. Specifically, the authors implicitly limited the regulatory region of an \"aging transcript\" by taking a region 1500bp upstream of the TSS. Naturally, this approach focuses on the promoter region but precludes the identification of many potentially relevant enhancer interactions. I wondered whether the authors have attempted to use epigenomic data, e.g., enhancer annotations for relevant blood cells, to detect such regulatory interactions.\n- I am not sure whether this is out of scope, but it would be interesting to learn more about the regulatory similarities of the identified aging genes. For instance, whether the set of multi-omic aging genes can be associated with a set of upstream regulators or whether there is an enrichment of transcription factor binding sites, e.g., by integrating TF-ChIP data.\n- The authors investigate the relationship of multi-omic aging genes and associated CpGs with mortality. In Figure 4c, Kaplan-Meier curves are shown for the top gene separately for the MGB-4K and GS data. I suggest the authors show the Kaplan-Meier curves for each gene/CpG in both data sets. For instance, adding the Kaplan-Meier plot for SATB1 (top gene in MGB-4K) obtained from the GS data and the Kaplan-Meier plot for TCF7 (top gene in GS) obtained from the MGB-4K data. Since many of the CpGs/genes identified in MGB-4K and GS are identical, I propose showing Kaplan-Meier curves for all of them, e.g., in a supplement.\nMinor recommendations:\n- I may have missed the detailed description of the PRC2-associated analyses in the Methods section. Otherwise, I suggest elaborating on this a bit more.\n- Although this information may be part of the code deposited on GitHub, I recommend explicitly stating which background model was used for GO analysis in the Methods section.", "reviewer_2": "The manuscript NCOMMS-24-55765-T entitled “Integrative epigenetics and transcriptomics identify aging genes in human blood” by Moqri et al. reports (i) the comparison, using six large cohorts, of the strongest correlations between DNAm levels and age with the strongest correlations between mRNA levels and age, replicating the previous observation that associations with the transcriptome are weaker; (ii) the observation that correlations between mRNA levels and age are weak for genes whose associated DNAm levels increase strongly with age, whereas correlations between DNAm and age are substantial for genes whose mRNA levels decrease strongly with age; (iii) evidence that genes showing an age-related increase in DNAm and an age-related decrease in mRNA (“multi-omic aging genes”) show more replicability across cohorts and are enriched in T-cell functions; (iv) evidence that aging CpGs are associated with the risk of dying with time and (v) aging CpGs remain associated with age when controlling for cellular composition.\nThe manuscript is well written and some results are of interest for the field. However, many aspects of the analyses are unclear and lack statistical support and justification. Please find below a list of comments and recommendations that should improve the manuscript.\nMajor comments\n- Most statements in the manuscript lack statistical support. I invite the authors to systematically report p-values of statistical tests supporting these statements. Examples: “we found that genes with the greatest decrease in expression levels with aging consistently featured significantly higher DNAm correlation with age”; “we found considerably higher and more consistent correlations with age across the datasets”; “multi-omic aging genes to be particularly enriched for T cell-specific genes”.… Another salient example is “age-associated DNAm data are considerably more replicable across cohorts than transcriptomic data”: to support this claim, the authors should report the proportions of CpGs and genes that show a significant correlation with age both in the discovery and replication cohorts, and then test if this proportion is significantly higher for DNAm.\n- The authors state that “multi-omic aging genes are strongly associated with clinically-relevant aging-related outcomes such as disease and mortality risk”. Actually, only DNAm levels at CpGs associated with aging genes are. Can the authors show volcano plots for expression levels of aging genes, as well as Peters et al.’s 1,497 genes? For comparison, I also suggest that the authors report volcano plots for the CpGs the most strongly associated with age, but that are not nearby aging transcripts. Finally, the authors are right that cellular composition can mediate the effects of age on DNAm or mRNA levels. I thus suggest that Cox regression models are adjusted on cellular composition.\n- The authors observed that correlations between DNAm and age are substantial for genes whose mRNA levels decrease strongly with age. Accordingly, Peters et al. have shown that “1,248 of the 1,497 age-associated genes had ≥1 CpG site with a significant Sobel test”, implying that DNAm mediates most of the associations between age and gene expression at aging transcripts. This result should be acknowledged and discussed in length. Notably,\n- It was previously demonstrated that CpG sites mediating age effects on gene expression are “less likely to reside in CpG islands or in promoters and more likely to be located in enhancers and insulators” (Peters et al., Nat Commun 2015). However, the authors of the current study defined “aging regulatory regions” as the 1,500 bases flanking the transcription start site of age-associated transcripts. Therefore, they may have missed most of the CpG sites that play a role in the regulation of genes whose expression is associated with age.\n- I concur with the authors that aging genes could be markers of T cell differentiation, mirroring the age-related shift of T cells from naïve to memory states. The authors rightfully tested whether CpGs associated with aging genes are associated with age when controlling for the cell proportions estimated by Houseman et al.’s or Horvath’s deconvolution methods. However, none of these methods estimate the proportion of the cell types that are most likely to mediate age effects: memory T cells. I thus invite the authors to report associations of DNAm and RNA levels with age for aging genes controlling for the proportion of the 12 cell types estimated by EPIC IDOL-ext (Salas et al., Nat Commun 2022).\n- The authors have used DNAm and RNA-seq from large cohorts of patients with age-related diseases (atherosclerosis, Parkinson). Therefore, disease status could mediate some of the strongest effects of age on DNAm levels or gene expression, whereas “aging genes” are supposed to be associated directly with age. Can the authors verify if their aging genes change when controlling for disease status (using multiple regression models)? Of note, the authors expect that “aging genes” could be used to “develop meaningful predictors for relevant aging-related outcomes”, but if these genes are directly associated such aging-related outcomes, interpretation would very different.\n- When the authors write “To validate our novel multi-omic aging genes, we aimed to perform an independent external validation using a novel cohort”, the reader does not know what are multi-omic aging genes, how they are defined exactly, their numbers... The authors cannot validate genes that are not yet identified.\n- It is actually unclear whether and why the authors have used MGB only to define the multi-omic aging genes. It would make more sense to discover aging transcripts in the large RNA-seq cohorts (not only the 1,497 transcripts from Peters et al. 2015), discover aging CpGs and genes in the large DNam and RNA-seq cohorts and then validate these genes in MGB, whose sample size is more modest. Related to this, when are “validated aging genes” used in the study?\n- The authors are right when they state that transcriptomic data show “higher noise, and susceptibility to batch effects” (Meyer et al., Aging Cell 2021). In the context of this study, were the RNA-seq data sets corrected for batch effects?\n- Are non-“aging genes” with age-associated expression enriched for T cell-specific genes?\n- Outlier values in DNAm and RNA-seq data can result in inflated correlations with age, particularly when the age distribution has long tails (e.g., 500FG, RA). Can the authors confirm their list of 106 aging genes with Spearman correlations?\n- It would be very informative to identify further the features that characterize aging genes. Are these genes enriched in binding site for specific TFs?\nMinor comments\n- I suggest that the authors clarify what is the added value of identifying “aging genes” compared to identifying genes directly associated with age-related diseases.\n- When the authors state their approach identifies “true aging-associated genes”, can they clarify what they mean?\n- Could the authors speculate further on the causes of the less reproducible effects of age on gene expression across cohorts, compared to DNAm? The CD248 example suggests that age has a non-linear effect on its mRNA levels.\n- Fig. 2: It could be informative to plot, for all genes, the correlation with age of gene expression against the correlation for associated DNAm.\n- Please italicize gene names in the text and figures.\n- In Fig 2c, it would be more informative to show results for PPMI instead of MESA 2010.\n- Fig 3c, please replace “correlation” by “correlation”\n- Table 1: Please add information on the tissue in which DNAm and RNA-seq were measured, and the percentage of diseased patients.\n- Methods: Replace “Pearrson” by “Pearson”, and “18,85911 participants” by “18,859 participants”.\n- Methods: Please describe in details the quality control filters applied on the new DNAm data."}, "author_response": "Point-by-Point Responses\nReviewer #1:\nIn \"Integrative epigenetics and transcriptomics identify aging genes in human blood,\" Moqri et\nal. present an integrative analysis of age-associated DNA methylation and RNA expression\nsignatures. The authors can pinpoint roughly 100 genes robustly associated with aging by\nanalyzing DNA methylation and RNA expression data from various cohorts. Importantly,\nenrichment analyses suggest that the identified genes may provide a better starting point for\nstudying the molecular mechanisms of aging in the human blood. The study achieves this\ngoal by generating an independent methylation data set to validate findings from previous,\npartially publicly available studies. Naturally, the generation of an independent RNA\nexpression data set of the same cohort would have been an added benefit, although not\nstrictly required, in my opinion, to substantiate the conclusions and claims of this study.\nConceptually, the manuscript introduces the term \"multi-omic aging genes.\" While it may not\nbe entirely new to seek genes exhibiting both transcriptional and regulatory changes during\nprocesses over time, the manuscript makes an essential case for analyzing multiple aging-\nassociated features simultaneously. It demonstrates that this approach helps to identify\nbiologically meaningful \"aging genes.\"\nThe manuscript is well-written, and the presented data substantiates all conclusions and\ninterpretations. Nevertheless, I have some suggestions that the authors may find helpful in\nrevising their manuscript.\n1. The authors may like to discuss the advantages and limitations of DNAm arrays. For\ninstance, arrays can measure only a fraction of all genomic CpGs, thus substantially limiting\nthe detection of age-associated DNA methylation alterations. On the other hand, arrays may\nexhibit higher reproducibility, e.g., compared to RRBS or WGBS.\nResponse: Thank you for the suggestion. DNA methylation arrays indeed measure only a\nsubset of all genomic CpGs, which limits our ability to detect the full spectrum of age-\nassociated DNAm alterations. Current generation arrays such as the Illumina\nMethylationEPIC v2.0 cover approximately 935,000 CpG sites, representing less than 5% of\nall CpG sites in the human genome. This selective coverage creates a potential bias in the\nidentification of age-associated DNAm signatures.\nHowever, arrays offer several advantages that make them particularly valuable for epigenetic\naging studies. First, they demonstrate exceptional technical reproducibility with typical\ncorrelation coefficients exceeding 0.98 between technical replicates. Second, the\nstandardized protocols and established analysis pipelines significantly reduce technical\nvariation between laboratories. Third, the relative cost-efficiency of arrays has enabled large-\nscale population studies that have been instrumental in defining robust epigenetic clocks.\nIn contrast, whole-genome approaches such as WGBS (Whole Genome Bisulfite Sequencing)\nand RRBS (Reduced Representation Bisulfite Sequencing) offer more comprehensive\ncoverage but face challenges in reproducibility and data processing. WGBS provides the most\ncomprehensive assessment of methylation status but requires significant sequencing depth to\nachieve accuracy comparable to arrays, substantially increasing costs for large-scale studies.\nRRBS offers a middle ground by enriching for CpG-rich regions but suffers from higher\ntechnical variability compared to arrays.\nAs suggested by this reviewer, we have added a brief section to our discussion, as follows:\n“While the DNA methylation array technology applied measures only a subset of\ngenomic CpG (approximately 935,000 CpG sites, representing 5%) in contrast to\nsequencing approaches, such as whole-genome bisulfite sequencing (WGBS), this\ntechnology offers several advantages for epigenetic studies, including technical\nreproducibility, standardized protocols and established analysis pipelines that\nsignificantly reduce technical variation between laboratories, and cost-efficiency”\n2. The current manuscript implicitly goes beyond integrating DNAm and RNA expression\nalready. Specifically, the authors show the enrichment of age-associated CpGs in regions\naffected by PRC2 binding. An increasing number of studies show that DNAm changes\noccurring during aging or in age-associated diseases are often associated with PRC2-\nassociated silencing. Also, because PRC2 plays a role in bivalently marked genes that\nfrequently play a role in development, I'd be curious whether the subset of multi-omic aging\ngenes with PRC2 association enriches any more specific functions.\nResponse: We appreciate the reviewer's insightful observation regarding PRC2 binding. Our\nenrichment analysis provides strong support for this connection, confirming our previous\nfindings (Moqri et al. 2024, see figures below). Specifically, our ChEA transcription factor\nanalysis revealed significant enrichment of PRC2-associated factors (p < 0.05) among our\nmulti-omic aging genes. This aligns with the growing literature linking age-associated DNA\nmethylation changes to PRC2-mediated regulation. Furthermore, our GO Biological Process\nanalysis identified significant enrichment of developmental processes, consistent with the\nbivalent marking characteristic of PRC2 target genes. These findings strengthen the\nmechanistic link between aging-associated epigenetic changes and developmental regulatory\nnetworks.\n3. The integration of DNAm and RNA expression data is straightforward but somewhat\nconservative. Specifically, the authors implicitly limited the regulatory region of an \"aging\ntranscript\" by taking a region 1500bp upstream of the TSS. Naturally, this approach focuses\non the promoter region but precludes the identification of many potentially relevant enhancer\ninteractions. I wondered whether the authors have attempted to use epigenomic data, e.g.,\nenhancer annotations for relevant blood cells, to detect such regulatory interactions.\nResponse: We appreciate the reviewer's comment regarding CpG sites associated with age-\nrelated gene expression. We have now assessed this using a systematic comparison of age\ncorrelation across different regulatory regions in aging-associated genes.\nAs shown in our analysis (added to the Supplement), we found that CpG sites in promoter\nregions demonstrate dramatically stronger age correlations compared to other genomic\nregions (p = 1.53e-36). The difference between enhancer regions and promoters was not\nstatistically significant (p = 1.13e-01), while enhancers showed modestly stronger correlations\nthan other regions (p = 2.43e-03). Please see the response to Reviewer 2 question 4 for more\ndetails.\n4. I am not sure whether this is out of scope, but it would be interesting to learn more about\nthe regulatory similarities of the identified aging genes. For instance, whether the set of multi-\nomic aging genes can be associated with a set of upstream regulators or whether there is an\nenrichment of transcription factor binding sites, e.g., by integrating TF-ChIP data.\nResponse: We thank the reviewer for this excellent suggestion. We have conducted\ncomprehensive transcription factor enrichment analyses using integration of ChIP-seq data\nthrough the ENCODE and ChEA databases revealed significant enrichment for specific\ntranscription factor binding (15 significant TF associations, p < 0.05). This analysis identified\nkey upstream regulators coordinating the expression of our multi-omic aging genes, notably\nmany of which are associated with the PRC2 complex (EZH2, SUZ12) or\nimmune/haematopoietic cell development (e.g., RUNX).\nThese findings reveal a coordinated regulatory network underlying age-associated epigenetic\nand transcriptional changes, rather than isolated alterations. This comprehensive regulatory\nanalysis strengthens our understanding of the molecular mechanisms driving aging-\nassociated genomic changes and provides new insights into potential therapeutic targets for\nage-related conditions.\n5. The authors investigate the relationship of multi-omic aging genes and associated CpGs\nwith mortality. In Figure 4c, Kaplan-Meier curves are shown for the top gene separately for the\nMGB-4K and GS data. I suggest the authors show the Kaplan-Meier curves for each\ngene/CpG in both data sets. For instance, adding the Kaplan-Meier plot for SATB1 (top gene\nin MGB-4K) obtained from the GS data and the Kaplan-Meier plot for TCF7 (top gene in GS)\nobtained from the MGB-4K data. Since many of the CpGs/genes identified in MGB-4K and\nGS are identical, I propose showing Kaplan-Meier curves for all of them, e.g., in a\nsupplement.\nResponse: We thank the reviewer for this suggestion. We have extended our KM analyses to\ninclude all the top CpG sites and included the results in the supplement, and copied below.\nMinor recommendations:\n- I may have missed the detailed description of the PRC2-associated analyses in the Methods\nsection. Otherwise, I suggest elaborating on this a bit more.\nResponse: We thank the reviewer for this suggestion. We have updated the figure legend to\nrefer the readers to the Method section for more details on the PRC2-associated analysis and\nadded the following details under the Methods:\n“Following Moqri et al. 2024, PRC2 binding regions were identified using EZH2 binding\nregions based on ChIP-Seq data from ENCODE Consortium (ENCFF109KCQ) and\nCpGs within these regions are marked as PRC2-associated sites”\n- Although this information may be part of the code deposited on GitHub, I recommend\nexplicitly stating which background model was used for GO analysis in the Methods.\nResponse: We thank the reviewer for this suggestion. We have updated the text to include\nthe following details about the GO analyses under the Methods:\n“Functional enrichment of genes associated with aging across modalities was\nperformed using Fisher’s exact test and gene terms from GO Biological Process, using\nall available genes in the ontology with clusterProfiler package in R. Redundant\nfunctions with high semantic similarity were filtered out with simplify function”\nReviewer #2:\nThe manuscript NCOMMS-24-55765-T entitled “Integrative epigenetics and transcriptomics\nidentify aging genes in human blood” by Moqri et al. reports (i) the comparison, using six large\ncohorts, of the strongest correlations between DNAm levels and age with the strongest\ncorrelations between mRNA levels and age, replicating the previous observation that\nassociations with the transcriptome are weaker; (ii) the observation that correlations between\nmRNA levels and age are weak for genes whose associated DNAm levels increase strongly\nwith age, whereas correlations between DNAm and age are substantial for genes whose\nmRNA levels decrease strongly with age; (iii) evidence that genes showing an age-related\nincrease in DNAm and an age-related decrease in mRNA (“multi-omic aging genes”) show\nmore replicability across cohorts and are enriched in T-cell functions; (iv) evidence that aging\nCpGs are associated with the risk of dying with time and (v) aging CpGs remain associated\nwith age when controlling for cellular composition.\nThe manuscript is well written and some results are of interest for the field. However, many\naspects of the analyses are unclear and lack statistical support and justification. Please find\nbelow a list of comments and recommendations that should improve the manuscript.\nMajor comments\n1. Most statements in the manuscript lack statistical support. I invite the authors to\nsystematically report p-values of statistical tests supporting these statements. Examples: “we\nfound that genes with the greatest decrease in expression levels with aging consistently\nfeatured significantly higher DNAm correlation with age”; “we found considerably higher and\nmore consistent correlations with age across the datasets”; “multi-omic aging genes to be\nparticularly enriched for T cell-specific genes”.… Another salient example is “age-associated\nDNAm data are considerably more replicable across cohorts than transcriptomic data”: to\nsupport this claim, the authors should report the proportions of CpGs and genes that show a\nsignificant correlation with age both in the discovery and replication cohorts, and then test if\nthis proportion is significantly higher for DNAm.\nResponse: We thank the reviewer for their valuable feedback regarding the need for\nstatistical support for our claims. We have now conducted comprehensive statistical tests\nthroughout our analysis and have updated the figures to include the relevant p-values and\nsignificance indicators.\nWe have analyzed the replicability of DNAm and transcriptomic data across cohorts using\nformal statistical testing. Our analysis confirms that age-associated DNAm data are\nsignificantly more replicable as aging CpGs showing higher correlations with age across\ndatasets compared to aging transcripts.\nWe have also validated the claim that \"genes with the greatest decrease in expression levels\nwith aging consistently featured significantly higher DNAm correlation with age.\" Statistical\ntesting across all three cohorts (MESA1, MESA2, and PPMI) demonstrated highly significant\ndifferences (Wilcoxon rank-sum test, p=4.24e-21 for MESA1, p=6.62e-08 for MESA2, and\np=1.99e-10 for PPMI), as indicated by the p-values displayed directly in figure 3f.\nRegarding our claim about multi-omic aging genes being \"particularly enriched for T cell-\nspecific genes,\" we have conducted formal enrichment analysis as shown in the graph below.\nThese results demonstrate significant enrichment for T-cell specific functions, including T cell\nreceptor signaling pathway (adjusted p-value < 0.01) and T cell differentiation pathways\n(Th1/Th2 and Th17 cell differentiation, adjusted p-values < 0.05).\nThroughout our study, we have systematically applied appropriate statistical testing to support\nour claims. We have updated all figures with statistical tests to include p-value annotations.\nWe believe these statistical tests address the reviewer's concerns and provide robust support\nfor the key claims in our manuscript.\n2. The authors state that “multi-omic aging genes are strongly associated with clinically-\nrelevant aging-related outcomes such as disease and mortality risk”. Actually, only DNAm\nlevels at CpGs associated with aging genes are. Can the authors show volcano plots for\nexpression levels of aging genes, as well as Peters et al.’s 1,497 genes? For comparison, I\nalso suggest that the authors report volcano plots for the CpGs the most strongly associated\nwith age, but that are not nearby aging transcripts. Finally, the authors are right that cellular\ncomposition can mediate the effects of age on DNAm or mRNA levels. I thus suggest that\nCox regression models are adjusted on cellular composition.\nResponse: We appreciate the reviewer's comment regarding the association between multi-\nomic aging genes and aging-related outcomes. We would like to clarify our approach and\nprovide additional evidence supporting our conclusions.\nAs the reviewer correctly notes, it is the DNAm levels at CpGs associated with aging genes\nthat show associations with clinical outcomes. We do not have access to gene expression\ndata in the MGB and GS cohorts used for outcome analyses, which prevents direct\ncomparison of expression-based associations. Instead, we conducted an alternative analysis\nto address this question.\nTo address the reviewer's concern, we selected the top 1,000 CpGs from MGB500 that\nshowed positive age correlation (DNAm gain), without pre-filtering based on association with\naging-related genes. We then generated volcano plots for mortality prediction in both GS and\nMGB4K cohorts. All Cox regression models were adjusted for age and 12 different cell type\ncompositions to account for the cellular heterogeneity concerns raised by the reviewer.\nThe volcano plots demonstrate that CpGs selected based on our multi-omic approach show\nstronger and more consistent associations with mortality across both cohorts compared to\nCpGs selected based solely on age correlation. The effect sizes for the non-gene-informed\nCpGs were notably smaller, with fewer sites reaching statistical significance after multiple\ntesting correction.\nThis cross-cohort validation supports our conclusion that integrating transcriptomic\ninformation with epigenetic data identifies more clinically relevant CpG sites. The adjustment\nfor cellular composition in our models further ensures that these associations are not merely\nreflecting age-related changes in blood cell proportions.\nThese findings reinforce our central claim that the multi-omic approach yields more robust\nbiomarkers of biological aging with stronger links to health outcomes than approaches using\neither data type in isolation.\n3. The authors observed that correlations between DNAm and age are substantial for genes\nwhose mRNA levels decrease strongly with age. Accordingly, Peters et al. have shown that\n“1,248 of the 1,497 age-associated genes had ≥1 CpG site with a significant Sobel test”,\nimplying that DNAm mediates most of the associations between age and gene expression at\naging transcripts. This result should be acknowledged and discussed in length.\nResponse: Thank you. We have now inserted a paragraph in the discussion section:\n“Our observation that genes with strong age-related decreases in expression also tend\nto show substantial age-associated DNAm changes is consistent with some prior\nfindings. Peters et al. (Nat Commun 2015) suggested that DNAm mediates the majority\nof associations between age and gene expression, with 1,248 of 1,497 age-associated\ngenes showing at least one CpG site with a significant Sobel test. However, there are\nnotable differences between our analyses. While they reported that such mediating\nCpGs were more likely to reside in enhancers and less likely in promoters which is in\ncontrast with our own findings, our analysis moreover focused specifically on CpG sites\nshowing robust methylation gain with age, a process known to preferentially affect\npromoter regions. Additionally, our study sought to identify “multi-omic aging genes”\nwhere both expression and DNAm shift with age, rather than test mediation. Thus,\nwhile the mechanistic focus differs, our findings are complementary and highlight the\ncomplexity of DNAm regulation across genomic contexts during aging.”\n4. It was previously demonstrated that CpG sites mediating age effects on gene expression\nare “less likely to reside in CpG islands or in promoters and more likely to be located in\nenhancers and insulators” (Peters et al., Nat Commun 2015). However, the authors of the\ncurrent study defined “aging regulatory regions” as the 1,500 bases flanking the transcription\nstart site of age-associated transcripts. Therefore, they may have missed most of the CpG\nsites that play a role in the regulation of genes whose expression is associated with age.\nResponse: We appreciate the reviewer's comment regarding CpG sites associated with age-\nrelated gene expression. We have now assessed this point through a systematic comparison\nof age correlation across different regulatory regions in aging-associated genes.\nAs shown in our analysis (added to the supplement, and figure copied below), we found that\nCpG sites in promoter regions demonstrate dramatically stronger age correlations compared\nto other genomic regions (p = 1.53e-36). The difference between enhancer regions and\npromoters was not statistically significant (p = 1.13e-01), while enhancers showed modestly\nstronger correlations than other regions (p = 2.43e-03).\nThis finding complements but also somewhat differs from Peters et al. (Nat Commun 2015),\nwho reported that CpG sites mediating age effects were \"less likely to reside in CpG islands\nor in promoters and more likely to be located in enhancers and insulators.\" It is important to\nnote that our study's focus is fundamentally different. While Peters et al. identified CpG sites\nthat mediate the effect of age on gene expression through Sobel mediation tests, our study\nspecifically aimed to identify robust \"multi-omic aging genes\" by integrating epigenetic and\ntranscriptomic age-related changes.\nMoreover, we specifically focused on age-related methylation gain, which has been shown in\nnumerous studies to preferentially affect promoter regions. The strong statistical significance\nwe observed (p = 1.53e-36) validates our methodological choice to focus on the 1,500 bases\nflanking transcription start sites as \"aging regulatory regions.\"\nThis approach has yielded a set of multi-omic aging genes that demonstrate greater\nreplicability across diverse populations and stronger associations with aging-related outcomes\ncompared to genes identified using either epigenetic or transcriptomic data alone, as detailed\nin our manuscript.\n5. I concur with the authors that aging genes could be markers of T cell differentiation,\nmirroring the age-related shift of T cells from naïve to memory states. The authors rightfully\ntested whether CpGs associated with aging genes are associated with age when controlling\nfor the cell proportions estimated by Houseman et al.’s or Horvath’s deconvolution methods.\nHowever, none of these methods estimate the proportion of the cell types that are most likely\nto mediate age effects: memory T cells. I thus invite the authors to report associations of\nDNAm and RNA levels with age for aging genes controlling for the proportion of the 12 cell\ntypes estimated by EPIC IDOL-ext (Salas et al., Nat Commun 2022).\nResponse. Thank you for this excellent point. In collaboration with the developer of IDOL\n(who has now joined as a co-author of this work), we have included the results for 12 cell\ntypes in the revision which is in line with our findings with 6 cell types. We have included\nthese figures in the supplement and updated the text accordingly as follows:\n“Significant associations of methylation levels at these CpGs with age remained after\nthe adjustment (Fig. 5e), indicating that age-related changes in these CpGs are\nindependent of changes in cell composition. The results are consistent when adjusting\nfor cell type compositions using 12 cell-type deconvolution (see Supplement)”\n6. The authors have used DNAm and RNA-seq from large cohorts of patients with age-related\ndiseases (atherosclerosis, Parkinson). Therefore, disease status could mediate some of the\nstrongest effects of age on DNAm levels or gene expression, whereas “aging genes” are\nsupposed to be associated directly with age. Can the authors verify if their aging genes\nchange when controlling for disease status (using multiple regression models)? Of note, the\nauthors expect that “aging genes” could be used to “develop meaningful predictors for\nrelevant aging-related outcomes”, but if these genes are directly associated such aging-\nrelated outcomes, interpretation would very different.\nResponse: This is an important point and we fully agree with this reviewer that disease states\ncan modulate epigenetic features.\nAll of our integrative aging genes are selected from the initial set of transcriptomic aging\ngenes identified using meta analyses of large diverse generally healthy or community-based\ncohorts by Peters et al. 2015. In addition, we have selected our epigenetic datasets from large\ncohorts of mostly healthy individuals and excluded diseased subjects from these cohorts. In\nparticular our reference DNA methylation dataset, MGB500, is selected from healthy controls\ndonors in the MGB Biobank, MESA participants were recruited from a generally healthy\npopulation (without any cardiovascular diseases), and Generation Scotland participants are\ngenerally healthier than the population (Navrady et al. 2017). We have now also included\nadditional details about the cohorts (in addition to the existing description of MGB4k and GS)\nand clarified the exclusion of diseased samples from specific cohorts, including the following:\n“MGB500 samples were selected from healthy control donors of diverse ages from the\nMGB Biobank. The donors were recruited from a major metropolitan academic medical\ncenter, were roughly balanced between male and female, and were generally\nrepresentative of the racial/ethnic distribution of the local area (Fig. 3a).\nMESA cohort is a population-based sample of an initial 6,814 men and women aged\n45-84 without known cardiovascular disease. Thirty-eight percent of the recruited\nparticipants were White, 28 percent African American, 22 percent Hispanic, and 12\npercent Asian, predominantly of Chinese descent. Participants were recruited from six\nfield centers across the United States: Wake Forest University, Columbia University,\nJohns Hopkins University, University of Minnesota, Northwestern University and\nUniversity of California - Los Angeles. (mesa-nhlbi.org).\nPublic DNAm and RNA-seq data were selected among commonly used datasets used\nin epigenetic or transcriptomic studies of aging from the generally healthy populations,\nwith the exception of PPMI which was included to increase our multi-omic data.\nDiseased participants from the GC6 cohort were excluded from the analyses.”\n7. When the authors write “To validate our novel multi-omic aging genes, we aimed to perform\nan independent external validation using a novel cohort”, the reader does not know what are\nmulti-omic aging genes, how they are defined exactly, their numbers... The authors cannot\nvalidate genes that are not yet identified.\nResponse: This is an important clarification. Thank you for noting. We have updated the text\nto elaborate this as follows.\n“To demonstrate that our integrative analyses and the identified epigenetic-\ntranscriptomic aging genes are replicable, we performed an independent external\nvalidation using a novel cohort we specifically designed for omic analysis of aging\nbiomarkers across lifespan. We generated DNAm profiles for 500 individuals of diverse\nages from the Mass General Brigham (MGB) Biobank using the Illumina Infinium\nMethylationEPIC v2.0 array, which covers over 935,000 CpG sites enriched for\nregulatory regions and has been shown to exhibit high reproducibility.”\n8. It is actually unclear whether and why the authors have used MGB only to define the multi-\nomic aging genes. It would make more sense to discover aging transcripts in the large RNA-\nseq cohorts (not only the 1,497 transcripts from Peters et al. 2015), discover aging CpGs and\ngenes in the large DNam and RNA-seq cohorts and then validate these genes in MGB, whose\nsample size is more modest. Related to this, when are “validated aging genes” used in the\nstudy?\nResponse: The MGB dataset was utilized as our primary discovery cohort despite its modest\nsize because it uniquely offers Epic2 array data with exceptional age range coverage and\nmatched multi-omic profiles from the same individuals. This platform provides superior\nregulatory element coverage and expanded CpG detection critical for comprehensive\nepigenetic aging analysis. We acknowledge the methodological considerations regarding\ncohort size and have added a statement to the discussion:\n“We acknowledge the MGB dataset's modest sample size as a limitation, though its\nEPIC2 platform provided superior CpG coverage critical for comprehensive epigenetic aging\nanalysis. Future studies should integrate EPIC2 technology with larger cohorts to validate our\nfindings and potentially identify additional aging biomarkers that may have been missed due\nto power constraints.”\n9. The authors are right when they state that transcriptomic data show “higher noise, and\nsusceptibility to batch effects” (Meyer et al., Aging Cell 2021). In the context of this study,\nwere the RNA-seq data sets corrected for batch effects?\nResponse: This is an important point and we indeed normalize expression data in Figure 1(c)\nwhere we analyse and combine multiple datasets. We have added a statement to the\nnormalization to the legend of Figure 1:\n“Gene expression data was normalized using median normalization”\nThe remaining analyses are done on individual RNA-seq data separately, hence batch effects\nacross datasets become less relevant. We have also adjusted the text in the methods as\nfollows:\n“We obtained previous RNA-seq datasets as raw gene expression counts. Count data\nwere normalized using log transformation on (raw counts + 1) and median adjusted to\naccount for batch effects when merging datasets. The majority of analyses was\nconducted within individual datasets rather than merging multiple datasets.”\n10. Are non-“aging genes” with age-associated expression enriched for T cell-specific genes?\nResponse: Thank you for this thoughtful question. We performed enrichment analyses for all\nthe genes with age-associated expression identified by Peters’ et al and found that they are\nnot in particular enriched for T cell-specific genes as shown below.\n11. Outlier values in DNAm and RNA-seq data can result in inflated correlations with age,\nparticularly when the age distribution has long tails (e.g., 500FG, RA). Can the authors\nconfirm their list of 106 aging genes with Spearman correlations?\nWe thank the reviewer for this methodological point about outlier effects on correlations. Our\nanalysis addressed this concern by calculating both Pearson and Spearman correlations for\nall CpGs with age in the MGB500 cohort. The high concordance between methods\n(correlation = 0.970) demonstrates that our findings are robust regardless of correlation\nmethod used. All 106 aging genes in our final list show significant correlations using both\nmethods (p < 0.05), confirming the reliability of our results independent of potential outlier\neffects.\n12. It would be very informative to identify further the features that characterize aging genes.\nAre these genes enriched in binding site for specific TFs?\nWe thank the reviewer for this excellent suggestion. We have conducted comprehensive\ntranscription factor enrichment analyses. Integration of ChIP-seq data through the ENCODE\nand ChEA databases revealed significant enrichment for specific transcription factor binding\n(15 significant TF associations, p < 0.05). This analysis identified key upstream regulators\ncoordinating the expression of our multi-omic aging genes, notably many of which are\nassociated with the PRC2 complex (EZH2, SUZ12) or immune/haematopoietic cell\ndevelopment (e.g., RUNX).\nThese findings reveal a coordinated regulatory network underlying age-associated epigenetic\nand transcriptional changes, rather than isolated alterations. This comprehensive regulatory\nanalysis strengthens our understanding of the molecular mechanisms driving aging-\nassociated genomic changes and provides new insights into potential therapeutic targets for\nage-related conditions.\nMinor comments\n- I suggest that the authors clarify what is the added value of identifying “aging genes”\ncompared to identifying genes directly associated with age-related diseases.\nResponse. Thank you for this excellent suggestion. We have included this in the discussion\nas follows.\n“We focus on “aging genes” rather than genes associated with specific age-related\ndiseases (e.g., Alzheimer’s disease, cardiovascular disease) to capture upstream,\nsystemic biological processes that drive disease susceptibility across multiple diseases\nwith age and may underlie a shared vulnerability to a spectrum of age-related\nconditions rather than disease-centric processes.”\n- When the authors state their approach identifies “true aging-associated genes”, can they\nclarify what they mean?\nResponse: As we demonstrate, there are many DNA methylation changes without a\nconcordant alteration in gene expression. Here, we identify genes that are jointly altered with\nboth gene expression and DNA methylation, suggesting a more relevant functional role. We\nhave clarified this in the discussion, as follows:\n“These analyses demonstrate that integration of multiple omics modalities has the\nability to refine the identification and validation of true aging-associated genes, where\nboth DNA methylation and gene expression are jointly altered, pointing towards\nfunctional shifts”\n- Could the authors speculate further on the causes of the less reproducible effects of age on\ngene expression across cohorts, compared to DNAm? The CD248 example suggests that\nage has a non-linear effect on its mRNA levels.\nResponse: We speculate that the less reliable reproducibility of the association between age\nand gene expression compared to DNA methylation is partially due to the challenges with\nmeasuring transcription levels compared to DNA methylation states. For example many RNA\nmolecules have relatively short life-time while most DNA methylation states are considered to\nbe relatively more stable.\n- Fig. 2: It could be informative to plot, for all genes, the correlation with age of gene\nexpression against the correlation for associated DNAm.\nResponse. Thank you for this suggestion. We have checked the relationship between the\nage correlation of gene expression and the DNA methylation and overall, there does not\nappear to be a strong association between these two correlations as shown in the figure\nbelow.\n- Please italicize gene names in the text and figures.\nResponse. Thank you for noting this. We have checked all the gene and transcript names\nand Italicized the names in the text when referring to genes. We have also taken a note of this\nand will make sure to Italicize all the gene names in the figure vectors that we will produce for\nthe publication.\n- In Fig 2c, it would be more informative to show results for PPMI instead of MESA 2010.\nResponse. In Fig 2c, the results in fact show data for all three cohorts (MESA 2000,\nMESA 2010, and PPMI).\n- Fig 3c, please replace “correlation” by “correlation”\nResponse. Thank you for noting this. We have corrected the typo in Fig 3c.\n- Table 1: Please add information on the tissue in which DNAm and RNA-seq were measured,\nand the percentage of diseased patients.\nResponse. Thank you for your suggestion. As mentioned in the text, all the samples are from\nblood and blood sample types are specified in Table 1, under “cell type”. We have now also\nadded additional details about the participants in each cohort, including the details stated\nunder point #6.\n- Methods: Replace “Pearrson” by “Pearson”, and “18,85911 participants” by “18,859\nparticipants”.\nResponse. Thank you for noting this. We have corrected these typos in the Methods.\n- Methods: Please describe in details the quality control filters applied on the new DNAm data.\nResponse. Thank you for this excellent point. We have added the following details under the\nMethod.\n“.IDAT files were preprocessed using standard parameters in the R SeSAMe package,\nversion 1.22.1, using the recommended processing for EPICv2 data, with the “QCDPB”\nparameter (qualityMask, inferInfiniumIChannel, dyeBiasNL, pOOBAH, noob). Samples\nand CpGs with more than 10 percent missing values were excluded from the\ngenerated beta value matrix.”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed all my comments and suggestions in full.", "reviewer_2": "The authors have satisfactorily addressed most of my comments.\n- My only remaining concern is the presence of batch effects within cohorts. Can the authors verify that first PCs of the RNA-seq and DNAm data do not reflect plate effects, beadchip effects or variation in sequencing coverage?\n- L. 123: Please italicize gene names.\n- Should the authors remove sentence in L. 130-131?\n- Resolution of Suppl. Fig. 12 should be improved. Could the authors clarify which CpGs were selected for analysis shown in panel b?\n- The Discussion is reported twice."}, "author_response": "Point-To-Point Responses to Reviewers' Comments\nReviewer #1 (Remarks to the Author):\nThe authors have addressed all my comments and suggestions in full.\nWe thank the Reviewer #1 for reviewing our manuscript.\nReviewer #2 (Remarks to the Author):\nThe authors have satisfactorily addressed most of my comments.\nWe thank the Reviewer #1 for reviewing our manuscript.\n- My only remaining concern is the presence of batch effects within cohorts. Can the authors\nverify that first PCs of the RNA-seq and DNAm data do not reflect plate effects, beadchip\neffects or variation in sequencing coverage?\nWe have performed additional analysis for the quality control in our MGB cohort and observed\nthat the plate numbers are not separated by the PC1.\nDNA methylation data from six independent cohorts (MESA, PPMI, Grady, MGB, GenoAge,\nand RA) were analyzed to investigate the distribution patterns of Type I and Type II CpG sites\ncomparing to PC1. The Illumina Infinium MethylationEPIC array utilizes two distinct probe\ndesign types for methylation detection. Type I probes employ two separate probes (one for\nmethylated and one for unmethylated cytosines) that bind to different sequences after bisulfite\nconversion, providing higher precision but covering fewer sites. Type II probes use a single\nprobe design that interrogates both methylated and unmethylated states at the same locus,\noffering broader genome coverage but with potentially different technical characteristics. The\nIllumina Infinium MethylationEPIC array manifest file (v1.0 B5) was used to classify CpG sites\nbased on their probe design type. Of the 865,918 CpG sites on the array, 862,927 sites had\ntype information available, comprising 142,137 Type I probes (16.5%) and 720,790 Type II\nprobes (83.5%). We show that there is no strong bias of Type 1 and Type 2 probes based on\nthe first PC loading.\nWe have added this part to our Methods section:\nQuality control analysis examined the distribution of the plate number as well as Type I and\nType II probe loadings on the first principal component across all cohorts to verify absence of\nsystematic technical bias (Supplementary Fig. 13).\nAll RNA-seq datasets used in this study were already quality-controlled and batch-corrected\nby the original data producers before we accessed. As we did not have access to the raw\nsequencing data to independently evaluate technical variation in sequencing coverage, we\nrelied on the preprocessing pipelines implemented by the respective consortiums.\nTOPMed (MESA) datasets underwent standardized laboratory methods with centralized\nprocessing using single pipelines to minimize batch, and PPMI RNA sequencing data were\nuniformly processed with standardized pipelines from RNA isolation through gene detection to\ncontrol technical. The original publications for these datasets report comprehensive batch\neffect correction procedures, including empirical Bayes methods and variance stabilization\napproaches appropriate for RNA-seq count data.\nOur analysis approach of conducting the majority of analyses within individual cohorts rather\nthan merging datasets further minimizes cross-cohort technical artifacts, with robust cross-\ncohort validation serving as an additional control against batch-driven findings.\n- L. 123: Please italicize gene names.\nWe have italicized gene names.\n- Should the authors remove sentence in L. 130-131?\nWe have removed it.\n- Resolution of Suppl. Fig. 12 should be improved. Could the authors clarify which CpGs were\nselected for analysis shown in panel b?\nWe have updated Supplementary Figure 12 to improve the resolution. The graph includes all\nthe CpGs we identified in our analysis (around the TSS of the 104 aging genes). We've also\nadded this description to the legend.\nSupplementary Fig. 12. Analysis of cell type composition and its impact on age-CpG\nassociations. a. Relationship between chronological age and estimated cell proportions for 12\ndifferent immune cell types. b. Scatter plot showing the association between age and CpG\nmethylation after correction for 12 cell types. Significant associations remained after adjustment,\nindicating that age-related changes in these CpGs are largely independent of changes in cell\ncomposition. All the CpG sites around the TSS of the 104 aging genes we previously identified are\nincluded in this analysis.\n- The Discussion is reported twice.\nWe have removed the extra copy of the discussion."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have satisfactorily addressed my remaining comments."}, "author_response": ""}]

Decision letter