Multi-omics network analysis reveals distinct stages in the human aging progression in epidermal tissue

Identification of latent age-associated molecular stages

To identify distinct stages in aging skin tissue, we examined 86 female subjects between 21 and 76 years. Subjects were chosen so that all ages were represented evenly and were required to be in good health. From each subject we sampled epidermal tissue from the subject’s volar forearms via the suction blister method. 31 of the original subjects were further re-invited for a longitudinal second measurement, which took place three years later (Figure 1). From the epidermal samples we generated gene expression and DNA methylation data, which were computationally integrated using the similarity network fusion approach. The resulting network, incorporating information from both transcriptomes and methylomes, was then used to identify hidden subtypes via unsupervised spectral clustering. The clustering revealed four distinct subgroups in our data with roughly equal sizes of 22 (cluster 1), 20 (cluster 2), 18 (cluster 3) and 26 (cluster 4) subjects, that captured the multi-omics similarity structure between the samples more clearly than either chronological or DNA methylation (DNAm) age (Figure 2A-2C). Association analysis to subject metadata showed that the clusters were significantly associated with chronological age (p = 5.8e-12, Figure 3A), whilst not being confounded by BMI (p = 0.71, Supplementary Figure 1A).

Figure 1. Study and analysis setup. Workflow diagram depicting the two-stage longitudinal study setup and the main steps of multi-omics data generation, integration and analysis.

Figure 2. Multi-omics similarity between subjects in the integrated network. (A) Heatmap visualization showing similarities between subjects in the fused multi-omics similarity network generated from gene expression and methylation data, with subjects ordered by increasing chronological age. (B) Same heatmap visualization of multi-omics similarity as in (A), with subjects ordered by increasing DNAm age. (C) Same heatmap visualization of multi-omics similarity as in (A) and (B), with the subjects ordered by the identified aging phases.

Figure 3. Biological age validation of the identified phases. (A) Boxplot showing chronological age distributions among the four identified aging phases. (B) Chronological age outliers among the aging phases, denoted as “old-like” for subjects that appeared to prematurely cluster into a higher aging phase, and “young-like” for subjects that were classified into a lower aging phase relative to their chronological age. (C) Boxplot showing the deviation of DNAm from chronological age based on aging phase outlier status, revealing a divergence in DNAm aging rate for aging phase outliers. Statistical significance determined using pairwise T-tests. (D) Hallmark of Aging signal strengths in gene expression data, comparing chronological age groups to the biological aging phases. Shown are the adjusted p-values from Anova comparisons, testing the segregation of the groupings among gene set enrichment scores. Figure adapted from the original Hallmark of Aging publication [15]. (E) Longitudinal validation after three-year period. The chord diagram shows aging phase classification of re-invited subjects at both time points, with phase transitions highlighted in red.

Latent stages associate more strongly with DNAm age rather than chronological age

As the unsupervised and purely data-driven clustering had identified groupings with strong association to chronological age, we explored the possibility that the clusters might capture hidden stages in the aging progression. We hypothesized that if this were so, the groupings ought to be more strongly associated with the actual biological aging state of our subjects, rather than their chronological age. To test this, DNAm ages of all subjects were calculated as previously described [19], to serve as a proxy measure for biological age (Supplementary Figure 1B). The comparison revealed a stronger association of the identified stages to DNAm age (p = 3.9e-13) as opposed to chronological age (p = 5.2e-12), strengthening the hypothesis that the clusters captured multi-omics aging stages.

Aging phase outliers are also biological age outliers in the sense of DNAm age

As subjects within the phases still presented considerable variation in chronological age, and the most proven approximation to biological age available to date is through the use of DNAm age, we further explored if the chronological outliers in the different aging phases were also outliers in the sense of DNAm age. We defined individuals as outliers for every aging phase (Figure 3B) if their chronological age either exceeded the 3^rd quartile by at least one third of the interquartile range (“young-like”) or subceeded the 1^st quartile by said amount (“old-like”). Analysis of the deviation of chronological to DNAm ages revealed that phase outliers were indeed biologically significantly younger (mean/median = -3.8/-4.5 y) or older (mean/median = +4.2/+4.3 y) than average and also than each other respectively (Figure 3C). Association testing using age-adjusted logistic regression models further revealed that subjects assigned to the “old-like” group were also significantly more likely to have reported frequent sun bathing in the questionnaires (p = 0.0336), delivering evidence of photoaging factoring into aging phase assignment. We further repeated the analysis using age estimates from a transcriptomic clock [19], which on average showed lower accuracy than its DNAm counterpart (Supplementary Figure 1C), in concordance with previous reports [19]. Association testing of phase outlier status to transcriptomic age revealed the same trends observed with DNAm age, albeit in this case without reaching statistical significance (Supplementary Figure 1D). Notably, the correlation of the two biological age markers was lower than their respective correlations to chronological age (Supplementary Figure 1E), indicating that the clocks capture at least partly independent features of aging, again underlining the importance of multi-omics approaches for aging research.

Aging phases show improved detection of Hallmark of Aging signals

The Hallmarks of Aging (HoA) describe nine main biological motives and processes that are believed to be driving the aging progression. We hypothesized that if the identified aging phases captured stages of biological aging, this ought to be reflected in gene expression patterns related to the known aging cornerstone processes, as summarized by the HoA. We therefore generated lists of genes involved in each of the nine HoA, by selecting GO and Reactome pathways which captured the essence of the respective hallmarks and combining them to novel gene sets (Supplementary Figure 2A). We then used the sets to test if the aging phases allowed better detection of HoA-enrichment signals than chronological age groups. ANOVA analyses using gene set enrichment scores indeed showed stronger discrimination based on aging stage for all HoA gene sets (Figure 3D). This further extends the evidence of the phases capturing biological multi-omics age to the level of gene expression.

Longitudinal validation of aging phases over three-year period

To assess if the phases could be longitudinally reproduced, 31 subjects from the original cohort were re-invited three years later for a second measurement. To assess the aging phase of the new samples, a random forest classification model was built on both expression and methylation features from data of the original cohort. The classifier demonstrated high accuracy (AUC = 0.95) in discriminating between the four phases in repeated cross-validation on the original data (Supplementary Figure 2B) and was subsequently used to predict the aging stages of all subjects at the second time point. For most subjects the aging phase did not change within the 3-year-period, indicating high stability of the identified groupings (Figure 3E). This finding is not unexpected, considering the time span of only three years past since original sampling, relative to the much larger average phase windows with a standard deviation of between 8.2-11.5 chronological years. Nonetheless five subjects could be observed migrating from one aging phase to another, all of them transitioning naturally along the age gradient to the next phase (Figure 3E). Notably, four of these five subjects were previously classified as chronological outliers at the upper end of their age phase (Supplementary Figure 2C).

Data-driven ranking of the Hallmarks of Aging along the phases reveals distinct succession patterns

In order to identify the most important biological motives driving the aging phases, we resorted to the use of machine learning models. For this we implemented a method based on classifiers that learn to distinguish between aging phases, whilst taking advantage of biological pathway information in the training process. The workflow consisted of a stepwise training of classifiers only on subsets of genes annotated to pathways to predict aging phase from gene expression. By restricting the training to these genes, the cross-validated accuracies of these classifiers allow the assessment of how well a given gene set enables the differentiation between aging phases, thus resulting in a score for each pathway’s relevance. This score can intuitively be interpreted as a measure of how important or predictive a gene set is to the grouping of interest. In order to derive a data-driven ranking of the Hallmarks of Aging along the aging phases, we performed this pathway predictivity analysis using the aforementioned HoA gene sets, calculating 100 permutations for each pathway model. The predictivity scores revealed a clear patterning of the HoA along the four phases that allowed a grouping of the hallmarks using hierarchical clustering (Figure 4A). Strikingly, the hallmarks clustered almost in the exact constellations postulated in their original description [15], namely into primary hallmarks (genomic instability, telomere attrition, epigenetic alterations and originally loss of proteostasis), antagonistic or secondary hallmarks (cellular senescence, deregulated nutrient sensing and mitochondrial dysfunction) and integrative hallmarks (altered intercellular communication and stem cell exhaustion). Our analysis did however reveal a divergence in the classification of the proteostasis-hallmark, which clustered more strongly with the group of integrative hallmarks. Examination of the HoA predictivity patterns based on the newly generated classification revealed that the predictivity peaks for the respective hallmark classes extracted through our analysis (Figure 4B) also precisely match the temporal manifestation sequence postulated in the original description of the hallmarks as well [15]: Namely, primary hallmarks peaked in aging phases 2 and 3, followed by a sharp drop in predictivity thereafter. Meanwhile the importance of the secondary hallmarks increased notably in aging phase 3. The integrative hallmarks, postulated to be emerging as a consequence of primary and secondary manifestation, increased slowly along the phases, while peaking in the late aging phases 3 and 4, again in concordance with the original postulation [15]. To our knowledge this is the first data-driven validation of the overarching sequence in which these cornerstones of aging manifest in human tissue.

Figure 4. Characterization of Hallmark of Aging predictivity within the aging phases. (A) Hierarchical clustering of the nine Hallmarks of Aging based on their gene set predictivity analysis along the four aging phases. Predictivity was determined using cross-validated random forest classifiers, trained to distinguish each of the aging phases from the others. (B) Predictivity of the Hallmark of Aging gene sets along the four aging phases, grouped into primary, secondary and integrative hallmarks. Statistical testing was performed using one-sided Wilcoxon tests. All predictivity scores were derived from 100 permutations.

Pathway predictivity analysis reveals a distinctly non-linear loss in pathway predictivity in old-age

As our analyses on the succession of the HoA already indicated a distinct shift in predictivities towards phase 4 of the identified aging stages, and in light of the recent publication of a sharp loss of signature identity for longevity-associated mTOR-signaling in human tissue around 60 years [14], we decided to expand our pathway analysis and explore the mid-to-late-life transition in more detail. For this we utilized the hallmark gene sets defined by the Broad Institute into the analysis, a set of conserved and highly refined gene sets created to improve pathway inference by reducing variance and gene overlap [20]. The predictivity analysis using these gene sets revealed a number of pathways that showed distinctly non-linear predictivity patterns along the four aging phases (Figure 5A). Notably, the most prominent global pattern was a sharp loss in predictivity observed in the transition from aging phase 3 to aging phase 4 in many pathways. In line with the loss of predictivity in nutrient sensing signaling hallmark observed in the HoA analysis and recent reports [14], mTOR-related signaling was among the pathways undergoing this distinct transition in the transition from phase 3 to phase 4 (Figure 5B), which, biological age outliers aside, matches the chronological age threshold of 60 identified in recent reports [14]. Other pathways exhibiting this pattern included oxidative phosphorylation and fatty acid metabolism, and notably also DNA repair pathways (Figure 5B). Exceptions from this trend included interferon and interleukin signaling, which increased steadily in predictivity along the phases, in line with the inflammaging theory of aging [21, 22], and the previously observed patterns in the HoA analysis. Apart from these exceptions, statistical analysis of all pathway predictivity signals between aging phases 3 and 4 still revealed a significant decrease in pathway predictivity, that is replicated using gene set enrichment analysis, also showing a distinct loss in pathway enrichment in transition to phase 4 (Figure 5C). As this finding potentially points to an increase in transcriptional noise, we investigated whether there was a change in the transcriptional similarity between subjects with transition into aging phase 4. For this we calculated pairwise correlations between the full transcriptomes of all subjects. In line with the results from the pathway analysis, we observed a significant drop in transcriptional similarity in the transition from phase 3 to phase 4. Notably, a similar effect can be observed in the methylation data, where a concomitant decrease in correlation between methylation profiles is observed (Figure 5D).

Figure 5. Global loss in pathway predictivity in the transition from mid- to late-life. (A) Heatmap showing the changes in pathway predictivity along the identified aging phases. The predictivities shown are the average predictivities calculated from 100 permutations for every pathway. (B) Scatterplots visualizing the changes in predictivity along the aging phases for selected pathways, several of which show distinctly non-linear patterns. (C) Overall loss in pathway predictivity observed in the transition from aging phase 3 to phase 4 is also detectable using gene set enrichment analysis. (D) Pairwise Pearson correlation between all subjects based on transcriptional and DNA methylation patterns.

Together these findings indicate a distinctly non-linear increase in biological noise in the transition from mid-to-late-life, likely to contribute to the deterioration of human tissue function in old age.

Recruiting

The study was performed in agreement with the recommendations of the Declaration of Helsinki and all test subjects provided written, informed consent. Subjects were recruited in the age range of 20 to 80 years, with equal numbers of participants within each decade. Subjects were required to be female, in good health and belonging to phototypes II or III according to the Fitzpatrick scale [33], to limit non-age related variability in the data. Exclusion criteria included tattoos or scars in the test area, pigmentation disorders, pregnancy and medication such as anti-histamines or anti-inflammatory drugs within two weeks prior to study start. A detailed listing of exclusion criteria can be found in the Supplementary. Participants were further required to complete a self-assessment questionnaire on age, weight, height, smoker status, sun bathing habits, as well as food and drinking habits upon study start.

Tissue sample preparation

The suction blister method applied in this study has been approved by the Ethics Commission of the University of Freiburg (general approval Dec 8, 2008; Beiersdorf AG No. 28857). Three suction blisters of 7 mm diameter were taken from the volar forearms of all subjects as previously described [34].

Nucleic acid extraction

Tissue samples were suspended in the respective lysis buffers for DNA or RNA extraction and homogenized using an MM 301 bead mill (Retsch). DNA was then extracted using the QIAamp DNA Investigator Kit (Qiagen) according to manufacturer’s instructions. RNA was extracted using the RNeasy Fibrous Tissue Mini Kit (Qiagen) according to manufacturer’s instructions.

Transcriptome sequencing

Transcriptome libraries were prepared using TruSeq Library Prep Kit (Illumina) and sequencing performed at 1x50 bp on Illumina’s HiSeq system to a final sequencing depth of 100 million reads per sample. Sequencing data was processed using a custom pipeline including Fastqc v0.11.7 [35] for quality control, Trimmomatic v0.36 [36] for trimming and Salmon v0.8.1 [37] for mapping and read quantification.

Array based methylation profiling

Methylation profiling was performed using Illumina 450k (first time point) and EPIC (second time point) arrays. In order to ensure comparability of measurements, EPIC arrays were computationally reduced to include only probes present on the original 450k array using the minfi package [38] in R [39]. Methylation data was processed in minfi using the funnorm normalization method.

Similarity network fusion and clustering

Prior to integration, the gene expression (log2 transformed transcripts per million) and CpG methylation data (M values) were batch corrected using the Combat algorithm [40] implemented in the sva package [41], following a feature selection step via filtering by median absolute deviation, retaining 10 % of the most informative features. The data was then integrated as previously described [18] using parameter settings of k = 10 (number of neighbors), t = 20 (number of iterations) and alpha = 0.5 (hyperparameter). Clustering on the fused network was performed via spectral clustering as previously described [18]. Measures used for the selection of cluster numbers were the eigen-gap statistic and rotation cost as proposed in the original method description [18], as well as visual inspection using heatmaps.

Age clock analyses

Analyses of DNAm and transcriptomic age were performed as previously described [19]. DNAm age was calculated from M values, whereas transcriptomic age was predicted based on log2 transformed transcripts per million.

Hallmark of aging gene sets

The HoA gene sets were generated from GO [42] and Reactome [43] gene sets by manually selecting matching pathways assigned to the nine Hallmarks of Aging [15]. A detailed list of genes annotated to each hallmark is provided in the Supplementary Material in .gmt format.

Enrichment analyses

Enrichment analyses were performed using the PLAGE algorithm based on singular value decomposition as described in [44] and implemented in the GSVA [45] R package.

Classification model to predict aging phase in longitudinal validation

To predict aging phase of re-invited subjects at the second time point, a random forest classifier was trained on the samples from the original cohort. Features were selected as the top 50 hits derived from differential gene expression analysis using DESeq2 [46] and differential methylation analysis using limma [47] from pairwise aging phase comparisons. The model was trained within the machine learning framework mlr [48], using the algorithm implemented in the original randomForest package [49]. Adjusted model hyperparameters included ntree = 1000 and $mtry\hspace{0.17em}=\hspace{0.17em}\sqrt{features}$. Accuracy of prediction was calculated as the area under the receiver operating characteristic curve (AUC) for multi-class comparisons, as implemented in the pROC package [50], and was derived from 5 x 5-fold repeated cross-validation.

Pathway predictivity analysis

Pathway predictivity was assessed using random forest pathway classifiers, constructed using the gene sets generated in this study and using the Hallmark Process [20] gene sets downloaded from the Molecular Signatures Database v6.2 [51]. The models were trained by restricting the molecular data to that of genes annotated within a given hallmark and trained to predict the aging phase of every sample. Predictivity was determined as the accuracy of correct classification derived from 5 x 5-fold repeated cross-validation for each pathway model, giving insight on how well genes within the gene set allow a discrimination between the phases and was thus used as a measure of importance of the respective pathway. Samples were stratified with respect to the target variable in the cross-validation process in order to avoid unbalanced proportions in any fold that might lead to bloated accuracy measures. Hyperparameters of all models were adjusted to ntree = 1000 and $mtry\hspace{0.17em}=\hspace{0.17em}\sqrt{numberofgenesinpathway}$. To determine the predictivity of the HoA stratified for each of four aging phases, the classifiers were separately trained in a one-against-all type of setup, learning to distinguish a phase from all the others. Modeling parameters and cross-validation were chosen as described above, and results for the four phases were aggregated afterwards.

General data analysis and visualization

Data analysis in R further included the usage of the package data.table [52], dplyr [53] and Hmisc [54] for data handling and general purpose functions, as well as the packages ggplot2 [55], ggpubr [56], ggsci [57], circlize [58] and pheatmap [59] for data visualization. Workflow diagrams were built using draw.io [60].