Development and ageing in multicellular organisms are highly intertwined processes. On the one hand, certain ageing-related phenotypes, such as presbyopia and osteoporosis (Luegmayr et al., 2004), are believed to represent the continuation of developmental processes into adulthood (Blagosklonny, 2006; de Magalhães and Church, 2005). Such cases of ‘runaway development’ or higher than optimal function during ageing (recognised as the hyperfunction theory of ageing; Gems and Partridge, 2013) may arise due to declined natural selection pressure failing to optimise expression regulation after sexual reproduction starts (Fisher, 1930; Medawar, 1953; Williams, 1957). Indeed, recent experimental studies in Caenorhabditis elegans show that senescence phenotypes promoted by insulin-IGF-1 signalling pathways support the hyperfunction theory (Lind et al., 2019; Ezcurra et al., 2018). On the other hand, molecular studies have also reported a reversal of the ageing transcriptome towards pre-adult levels in various contexts, including primate brain regions (Somel et al., 2010; Dönertaş et al., 2017; Colantuoni et al., 2011), and mouse liver and kidney (Anisimova et al., 2020). Studying the functional consequences of this reversal pattern in the ageing human brain, we previously interpreted it as an indication of loss of cellular identity in neurons, possibly exacerbated by a reduction in the relative frequencies of neurons (Dönertaş et al., 2017). Such changes, in turn, could be caused by the accumulation of stochastic damage at the genetic, epigenetic, and proteomic levels over an adult lifetime, causing deregulation of gene expression networks.

Several major questions remain. First, the prevalence of reversal phenotypes across tissues is unclear as most research has been conducted in the brain (Somel et al., 2010; Dönertaş et al., 2017). A second question pertains to the similarity of reversal-exhibiting genes and pathways across tissues. Ageing-related expression changes are partly shared among organs (Zahn et al., 2007), and reversal trends are also shared across different regions of the primate brain (Dönertaş et al., 2017). Distinct tissues might hence show parallel reversal patterns. Alternatively, as mammalian tissues diverge from each other during development in their transcriptome profiles (Cardoso-Moreira et al., 2019), one may hypothesise that during ageing tissues converge back towards similar transcriptome profiles. Such a putative late-age convergence phenomenon would be consistent with the notion of ageing-related cellular identity loss (Yang et al., 2019; Dönertaş et al., 2017). A final question concerns the mechanism behind the observed reversal trends at the bulk tissue level. Specifically, the contribution of cell-type composition and cell-autonomous changes to the reversals at the tissue level remains unexplored.

Documenting the reversal phenomenon is critical to better understand the proximate mechanisms of mammalian ageing, and its ultimate mechanisms, such as the stochastic disruption versus continued expression of developmental genes. However, such work has been limited by the scarcity of studies that include both development and ageing periods of the same organism and across different tissues. This work presents an age-series analysis of bulk transcriptome profiles of mice, including samples of four tissues across postnatal development and ageing periods covering the whole postnatal lifespan. Using this dataset, we study the prevalence, mechanisms, and functional consequences of the reversal phenomenon in different mouse tissues. We further test the related hypothesis of tissue convergence during ageing and investigate the contribution of cell-type composition and cell-autonomous changes.

We generated bulk RNA-seq data from 63 samples covering the cerebral cortex (which we refer to as cortex), liver, lung, and skeletal muscle (which we refer to as muscle) of 16 male C57BL/6J mice, aged between 2 and 904 days of postnatal age (Materials and methods). As mice reach sexual maturity by around 2 months (Tacutu et al., 2018), we treated samples from individuals aged between 2 and 61 days (n = 7) as the development series, and those aged between 93 and 904 days (which roughly correspond to 80-year-old humans; Flurkey et al., 2007) (n = 9) as the ageing series (Figure 1—figure supplement 1). The final dataset contained n = 15,063 protein-coding genes expressed in at least 25% of the 63 samples (one 904-day-old mouse lacked cortex data).

Tissues diverge during postnatal development

Consistent with earlier work (Brawand et al., 2011; Cardoso-Moreira et al., 2019), we found that variation in gene expression is largely explained by tissue differences, such that the first three principal components (PCs) separate samples according to tissue (ANOVA p<10^–20 for PC1–3, Figure 1—source data 1), with the cortex most distant from the others (Figure 1a). Meanwhile, PC4, which explains 8% of the total variance, displayed a shared age effect across tissues in development (Spearman’s correlation coefficient ρ = [-0.88, –0.99], nominal p<0.01 for each test; Figure 1b). Also, after the tissue effect was removed by standardisation, principal components analysis (PCA) showed a strong influence of age on the first two PCs, which explains 31% of the variance in total (Figure 1—figure supplement 2). We further observed higher similarity among tissues at the juvenile stage compared to the young-adult stage. In other words, distances between tissues increased with age (change in mean Euclidean distance among tissues with age during development in PC1–PC4 space ρ_dev = 0.99, p_dev = 1.5 × 10^–5, Figure 1—source data 1), which resonates with previous reports of inter-tissue transcriptome divergence during development (Cardoso-Moreira et al., 2019). This divergence pattern was also observed when PCA was performed with developmental samples only (days 2–61: change in mean Euclidean distance among tissues in PC1–PC4 space; ρ = 0.95, p=0.0008; Figure 1—figure supplement 3a and b).

Figure 1 with 15 supplements see all

Download asset Open asset

Figure 1—source data 1

Data summary, age-related expression patterns, and reversal patterns.

https://cdn.elifesciences.org/articles/68048/elife-68048-fig1-data1-v4.xlsx

Download elife-68048-fig1-data1-v4.xlsx

Inter-tissue divergence during development and convergence during ageing

We studied the inter-tissue divergence/convergence question using two approaches. In the first, we analysed how transcriptome-wide expression variation among tissues changes with age regardless of their age-related expression patterns in any particular tissue. To do this, for each individual, we calculated the coefficient of variation (CoV) across the four tissues for each commonly expressed gene (n = 15,063), which represents a measure of expression variation among tissues. Then, we assessed how such inter-tissue variation changes over the lifetime by calculating the Spearman’s correlation coefficient between CoV and age separately for development and ageing periods (correlation values for all genes are given in Figure 2—source data 1).

Using the CoV values calculated across all 15,063 genes (excluding one 904-day-old individual for which we lacked the cortex data), we observed a significant mean CoV increase in development (Spearman’s correlation coefficient ρ = 0.77, two-sided p=0.041), confirming that tissues diverge as development progresses (Figure 2a). Interestingly, during ageing, we observed a decrease in mean CoV with age, albeit not significant (ρ = −0.50, p=0.204, Figure 2a), suggesting that tissues may tend to converge during ageing. This was also supported by the PCA in which we observed a trend of ageing-associated decrease in mean Euclidean distance among tissues (using PC1–PC4 space with quantile-normalised data: ρ = −0.87, p=0.0026; with VST-normalised data ρ = −0.58, p=0.102, Figure 1—source data 1). We obtained the same divergence-convergence pattern by calculating the median CoV values for each individual instead of the mean (Figure 2—figure supplement 1). Figure 2b exemplifies this pattern of increasing and then decreasing CoV through lifetime for the gene displaying the strongest such signal.

Figure 2 with 20 supplements see all

Download asset Open asset

Figure 2—source data 1

All the data related to divergence-convergence (DiCo) pattern: age-related coefficient of variation (CoV) change of genes, pairwise tissue expression correlations, analysis of independent datasets; GSE34378 (Jonker et al.), GSE132040 (Schaum et al.), and GTEx.

https://cdn.elifesciences.org/articles/68048/elife-68048-fig2-data1-v4.xlsx

Download elife-68048-fig2-data1-v4.xlsx

We identified n = 9058 genes showing divergent trends among tissues in development based on their CoV change with age (without using a significance cutoff per gene). Among these, n = 4802 showed convergent trends in ageing, which we refer to as DiCo genes. We next studied the transition points between divergence and convergence by clustering genes showing the DiCo pattern (n = 4802) based on their CoV values (Figure 2—figure supplement 2). Notably, cluster 1, which shows a slightly delayed divergence starting after 8 days and peaks around 3 months, was associated with metabolic and respiration-related processes (FDR-corrected p-value<0.1), and cluster 5, which shows a relatively delayed convergence after 4 months, was enriched in categories related to vascular development (FDR-corrected p-value<0.1) (Supplementary file 4). To assess the contribution of different tissues to the DiCo pattern, we further clustered DiCo-displaying genes (n = 4802) based on their expression levels (Figure 2—figure supplement 3). Not surprisingly, the clusters with relatively higher expression levels of a tissue (e.g. muscle in cluster 9) were enriched in functional categories (FDR-corrected p-value<0.1) related to that tissue (e.g. muscle cell development) (Supplementary file 5).

We then studied DiCo at the single-gene level. We tested each gene for a significant CoV change in their expression levels (i.e. divergence or convergence) in development and ageing (Spearman’s correlation test with FDR-corrected p-value<0.1). We found that the ratio of divergent and convergent genes differed significantly between development (70% divergence among 2581 significant genes) and ageing (68% convergence among 62 significant genes) (Figure 2d and e). The same pattern was also observed without using significance cutoff (Figure 2—figure supplement 4). We also confirmed that this pattern is also observed with VST-normalised data (Materials and methods), and is thus not affected by the data preprocessing approach (Figure 2—figure supplement 14).

To our knowledge, inter-tissue convergence during ageing is a novel phenomenon. We first considered the possibility that convergence during ageing could be explained by heteroscedasticity which could arise due to increased inter-individual variability in gene expression during ageing (Somel et al., 2006). To test this hypothesis, we compared expression-age heteroscedasticity levels between two gene sets: (1) genes with the DiCo pattern and (2) genes showing divergent patterns throughout lifetime (divergent-divergent [DiDi, n = 4182]) for each tissue separately (Materials and methods). We did not observe any significant difference in heteroscedasticity between DiCo and DiDi genes in any of the tissues (two-sided Kolmogorov–Smirnov [KS] test, p>0.05 in all tissues, Figure 2—figure supplement 15), which suggests that heteroscedasticity due to increased inter-individual variability probably does not drive the observed age-related convergence during ageing. Visual inspection of gene expression clusters also suggested that the DiCo pattern is not particularly associated with nonlinear changes in gene expression with age (Figure 1—figure supplements 12–15).

In order to further verify the DiCo pattern, we used a second approach to test it in our mouse dataset. For each individual, we calculated correlations between pairs of tissues across their gene expression profiles. Under the DiCo pattern, we would expect pairwise correlations to decrease during development and increase during ageing. Among all pairwise comparisons, we observed a strong negative correlation during development (ρ = [-0.61, –0.9], nominal p<0.05 in five out of six tests), while during ageing, four out of six comparisons showed a moderate positive correlation (ρ = [0.16, 0.69], nominal p<0.05 in one out of six comparisons, Figure 2—figure supplement 5). Calculating the mean of pairwise correlations among tissues for each individual, we observed the same DiCo pattern (nominal p<0.05 for both periods, Figure 2—figure supplement 6).

The DiCo pattern indicates loss of tissue specificity during ageing

Potential explanations of the DiCo pattern involve two scenarios consistent with the age-related loss of identity: (1) decreased expression of tissue-specific genes in their native tissues or (2) non-specific expression of tissue-specific genes in other tissues. To test these predictions, we first identified tissue-specific gene sets based on relatively high expression of that gene in a particular tissue (cortex: 1175; lung: 839; liver: 986; muscle: 766 genes). We noted that tissue-specific genes show clear up-down reversal patterns, being mostly upregulated during development, and downregulated during ageing (Figure 3, 57–89%). The up-down reversal pattern was particularly strong among tissue-specific genes for the three of four tissues tested (OR = [1.65, 6.52], p<0.05 for each tissue except in liver: OR = 0.87, p=0.09, Figure 3—source data 1). Tissue-specific genes were also enriched among DiCo genes (Figure 3—source data 1, OR = 1.56, Fisher’s exact test p<10^–16).

Figure 3

Download asset Open asset

Figure 3—source data 1

Effect sizes for determination of tissue-specific genes, enrichment of divergence-convergence (DiCo), and reversal genes within tissue-specific genes.

https://cdn.elifesciences.org/articles/68048/elife-68048-fig3-data1-v4.xlsx

Download elife-68048-fig3-data1-v4.xlsx

We then tested our initial prediction that the DiCo pattern is related to tissue-specific genes losing their expression in their native tissue and/or gaining expression in non-native tissues during ageing. We first tested this hypothesis by considering all tissue-specific genes. We found a positive odds ratio between loss of expression in native tissue and gain in other tissues during ageing (OR = 5.50, Fisher’s exact test p=2.1 × 10^–129, Figure 4a). The same analysis conducted with only the DiCo genes yielded a much stronger association (OR = 74.81, Fisher’s exact test p=5.9 × 10^–203, Figure 4b). This suggests that loss of tissue-specific expression is observed across the transcriptome, with a particularly strong association among DiCo genes. Figure 4c–f exemplifies the expression trajectories of genes chosen from each group defined in Figure 4b.

Figure 4 with 2 supplements see all

Download asset Open asset

Figure 4—source data 1

Gene Set Enrichment Analysis (GSEA) result of divergence-convergence (DiCo) genes, DiCo enrichment with tissue-specific expression loss, age-related expression change correlations, and convergence overlaps among datasets.

https://cdn.elifesciences.org/articles/68048/elife-68048-fig4-data1-v4.xlsx

Download elife-68048-fig4-data1-v4.xlsx

We then asked whether genes displaying the DiCo pattern may be related to specific functional pathways or share specific regulators. Using GO, we searched for functional enrichment among convergent genes during ageing using developmentally divergent genes as the background (Materials and methods). We found enrichment for 184 GO Biological Process (BP) categories for the DiCo pattern (KS test, FDR-corrected p-value<0.1, Figure 4—source data 1) and summarised enriched categories by clustering them based on the number of genes they share. We then studied the trends of gene expression changes with age (without a significance cutoff) in each representative category for each tissue (Materials and methods) (Figure 4h; we provide detailed clustering for the categories in ‘Other GO’ Figure 4—figure supplement 1). On average, energy metabolism, mitochondria, and tissue function-related categories, as well as immune response-related categories, exhibit DiCo-type expression changes over time and across tissues, where temporal changes in different tissues occur in opposite directions. Notably, for the majority of representative GO categories, the lung had the most distinct expression patterns in both periods (Figure 4h, Figure 4—figure supplement 1).

Contrary to the functional enrichment results, we did not find any specific regulators (miRNA or transcription factors) associated with DiCo using the same background as above (at 235 tests for miRNA and 158 tests for TF, FDR-corrected p-value>0.1 for both tests) (Materials and methods), which suggests that DiCo pattern may not be driven by a limited number of specific regulators, but may instead be a transcriptome-wide phenomenon.

Changes in cellular composition and cell-autonomous expression can both explain the DiCo pattern

Ageing-related transcriptome changes observed using bulk tissue samples may be explained by temporal changes in cell-type proportions within tissues, by cell-autonomous expression changes, or both. To explore whether the observed inter-tissue DiCo patterns may be attributed to changes in cell-type proportions, we used published data from a mouse single-cell RNA-sequencing experiment (Tabula Muris Consortium, 2020). For each of the four tissues in our original experiment, we collected cell-type-specific expression profiles from 3-month-old young adult mice in the Tabula Muris Senis dataset. We deconvoluted bulk tissue expression profiles in our mouse dataset using the corresponding tissue’s cell-type-specific expression profiles by regression analysis (Materials and methods) and studied the relative contributions of each cell type to tissue transcriptomes and how these change with age. The analysis was performed with three gene sets; all genes (n = [12,492, 12,849]), DiCo (n = [4007, 4106]), and non-DiCo genes (n = [8485, 8743]). Studying these deconvolution patterns, we observed a weak but consistent trend involving the most common cell types in different tissues. For instance, analysing DiCo genes in the liver and lung, we found that the most common cell type’s contribution (hepatocyte in the liver, and bronchial smooth muscle cell in the lung) tends to increase during development (Spearman’s correlation coefficient ρ_liver = 0.95, ρ_lung = 0.81, nominal p<0.05). This contribution then decreases during ageing (ρ_liver = –0.77, ρ_lung = –0.86, nominal p<0.05) (Figure 5a, Figure 5—figure supplement 1). This pattern was also observed in muscle and cortex, albeit not significantly (Figure 5a, Figure 5—figure supplement 1). These changes most likely reflect shifts in cellular composition, some of which were demonstrated directly in mice using in situ RNA staining (Schaum et al., 2020). Repeating the analysis with non-DiCo genes resulted in highly similar patterns considering the most common cell types in tissues, except in muscle ageing in which the age-related decrease was significantly higher with DiCo genes than the non-DiCo genes (permutation test with resampling all genes, p_{skeletal-muscle-satellite-cell} = 0.04) (Figure 5a, Figure 5—figure supplement 1, Figure 5—figure supplements 2–5). These results indicate that the observed cellular composition changes may partly explain DiCo, although the influence of composition changes is not exclusive to genes displaying the DiCo pattern.

Figure 5 with 6 supplements see all

Download asset Open asset

Figure 5—source data 1

Cell-type proportion estimation and cell-autonomous changes using the Tabula Muris Senis dataset.

https://cdn.elifesciences.org/articles/68048/elife-68048-fig5-data1-v4.xlsx

Download elife-68048-fig5-data1-v4.xlsx

Next, we investigated the possible role of cell-autonomous changes in the DiCo pattern. Cell-autonomous changes could contribute to inter-tissue convergence during ageing in two ways. First, expression profiles of similar cell types shared across different tissues, such as immune cells, might converge with age. Another possible scenario, consistent with the notion of age-related cellular identity loss, is that the expression profiles of unrelated cell types, such as tissue-specific cell types in different tissues, converge with age. To test these scenarios, we first ordered the pairwise correlations between cell types in different tissues at the 3-month age group to determine the most similar and dissimilar cell types across tissues (Materials and methods). Then, we studied how these similarities (i.e. pairwise correlations) change with age (Figure 5b). Intriguingly, we found that pairs of similar cell types (i.e. those with the highest correlations) among tissues tend to become less similar with age (36/54 [67%] of pairwise comparisons, Figure 5—source data 1). On the contrary, the most distinct cell types (i.e. those with the lowest correlations) among tissues become more similar with age (45/54 [83%], Figure 5—source data 1). Repeating the analysis considering DiCo genes only yielded a similar trend (30/54 [56%] decrease in correlation among the most similar cell types, permutation test with resampling non-DiCo genes, p>0.1; and 47/54 [87%] increase in correlation among the most distinct cell types, permutation test, p>0.1). These trends are consistent with age-related cellular identity loss, and they suggest that cell-autonomous changes may also contribute to inter-tissue convergence during ageing, although further data and analyses would be needed to fully establish their validity.

Finally, we tested the possibility of intra-tissue convergence of cell types in the Tabula Muris Senis dataset by calculating expression variation among cell types using the CoV measure for each individual. However, we did not observe a consistent trend of increasing similarity among cell types within tissues from 3-month-old to 24-month-old mice (Figure 5—figure supplement 6).

Our findings confirm a number of ageing-associated phenomena identified earlier, while also revealing new patterns. First, we report parallel age-related expression changes among the four tissues studied, during development, as well as in ageing. The inter-tissue correlation distributions were modest and also comparable between development and ageing (Figure 1c). This last point may appear surprising at first glance, given the stochastic nature of ageing relative to development (Bahar et al., 2006; Martinez-Jimenez et al., 2017; Angelidis et al., 2019; Somel et al., 2006; Feser et al., 2010; Kim et al., 1996; Enge et al., 2017), and also given earlier observations that developmental expression changes tend to be evolutionarily conserved, while ageing-related changes much less so (Zahn et al., 2007; Somel et al., 2010). At the same time, when we consider that tissues diverge during development, and also that ageing is characterised by parallel expression changes among tissues related to damage response, inflammation, and reduced energy metabolism (Zahn et al., 2007; Yang et al., 2015), similar magnitudes of correlations during development and ageing may be expected.

Second, we verify the generality of the reversal pattern, that is, up-down or down-up expression change patterns across the lifetime, among distinct mouse tissues that include both highly mitotic (lung and liver) and less mitotic ones (skeletal muscle and cortex). Consistent with earlier observations in fewer tissues (Anisimova et al., 2020; Dönertaş et al., 2017), we find that about half the expressed genes display reversal in all cases studied. Importantly, expression reversal is not ubiquitous across all genes and our findings do not necessarily contradict the hyperfunction theory. Instead, we suggest that reversal is a common phenomenon that influences a notable fraction of the transcriptome and is a likely contributor to mammalian ageing.

Two observations here are notable. One is that reversal-displaying genes, especially those displaying the up-down pattern in each tissue, can be associated with tissue-specialisation-related pathways (e.g. morphogenesis) and tissue-specific functions (e.g. synaptic activity). The second observation is the lack of significant overlap among reversal genes among tissues. We thus hypothesised that reversals might be reflecting tissue specialisation during development (hence lack of overlap among tissues) and loss of specialisation during ageing. These processes could manifest themselves as inter-tissue divergence and convergence patterns over lifetime. We indeed observed that the up-down reversal pattern is enriched in tissue-specific genes, except in the liver. Studying inter-tissue similarity across mouse lifespan, we further found that the four tissues’ transcriptomes diverged during postnatal development, and we further detected a trend towards inter-tissue convergence during ageing. We then further investigated this phenomenon through different approaches: (1) by studying overall trends using PCA, (2) by analysing transcriptome-wide trends of inter-tissue CoV without considering gene-wise significance cutoffs, (3) by focusing on genes with significant age-related changes in inter-tissue CoV, (4) by studying age-related changes in pairwise tissue correlations, (5) by analysing different cell types using scRNA-seq data, and (6) by repeating the same analysis using independent mouse and human ageing datasets. The patterns we found were mostly consistent with inter-tissue convergence, but the majority of transcriptome-wide results were associated with low ESs, and some were not statistically significant. Importantly, all significant results suggested convergence during ageing. We therefore conclude that (1) developmental inter-tissue divergence does not continue into ageing and (2) convergence during ageing may be common although possibly not ubiquitous.

The weakness of the inter-tissue convergence signal per dataset and the limited overlap between convergent gene sets among datasets could have multiple reasons. These include the low signal-to-noise ratios characterising ageing-related expression patterns, the lack of old-age individuals in our mouse dataset (>3-year-old mice) and the GTEx dataset (>90-year-old humans), limited overlap of tissues between our mouse dataset (cortex, liver, lung, and muscle) and the Jonker et al. dataset (cortex, liver, lung, spleen, kidney), as well as differences in ageing patterns between species or between sexes. Further research involving larger sample sizes and diverse species is needed to confirm the generalisability of the observations.

Finally, we report a number of interesting observations on DiCo. We determine that tissue-specific genes tend to be downregulated in the tissues that they belong to during ageing, while non-tissue-specific genes are upregulated, which was confirmed by all external datasets (Figure 4c). Second, using deconvolution, we infer that cell types most common in a tissue (e.g. hepatocytes in the liver) tend to increase in frequency during development, but then decrease in frequency during ageing, as also shown recently using immunohistochemistry in a number of mouse tissues (Schaum et al., 2020). Accordingly, the DiCo phenomenon may at least partly be explained by shifts in cellular composition. This is intriguing as both highly mitotic and low mitotic tissues share this trend, indicating that an explanation based on stem cell exhaustion may not be applicable here. Third, we find increased expression similarity between distinct cell types in different tissues during ageing, but decreased similarity between similar cell types. Cell-autonomous expression changes, therefore, likely also contribute to the DiCo phenomenon. We note that higher expression variability among cells at old age (Hernando-Herraez et al., 2019; Enge et al., 2017) could also lead to inter-tissue convergence during ageing. A fourth interesting observation was the absence of significant enrichment for specific transcription factor or microRNA targets among DiCo genes. This result may not be surprising if inter-tissue convergence is mostly driven by stochastic damage accumulation, such as loss of epigenetic marks. It is also possible that instead of specific regulators their interaction and cooperativity are associated with the DiCo. Future experimental studies could test both mechanistic aspects and functional link to tissue specificity.

We also note two major limitations of our study. One is related to the fact that our dataset represents bulk tissue samples, which may suffer from infiltration of foreign cell types into tissues. Indeed, one of the external datasets, Schaum et al., included samples from perfused mice (Schaum et al., 2020), and we did not find support for the transcriptome-wide convergence during ageing, even though the association between tissue identity loss and convergence was also evident. The scRNA-seq dataset we analysed further suggested that DiCo is associated with tissue-specific genes and not immune- or blood-related categories, but we still cannot rule out possible infiltration artefacts that may affect our results. A second limitation is related to ageing being highly sex-dimorphic in mammals (Yuan et al., 2012; Sampathkumar et al., 2020). Hence, in-depth analysis of sex specificity of the DiCo pattern could be relevant. Our mouse dataset included only male mice, while that of Jonker et al. was female-only. The fact that both revealed DiCo patterns suggest DiCo is not particular to one sex, but there could still exist sex-specific effects. In fact, when we analysed DiCo among human male and female individuals in the GTEx dataset separately, we observed slightly stronger inter-tissue convergence among ageing females than in males, although the GTEx male samples have also a drastically narrower age range (Figure 2—figure supplement 16). Accordingly, the prevalence of DiCo among humans and sexes waits to be determined.

Despite the open questions that remain, our results consistently support a model where ageing mammals suffer from loss of specialisation at the tissue level, and possibly also at the cellular level, which are observed as expression reversals and the newly discovered DiCo phenomenon we report here.

Sequencing data generated for this study have been deposited in GEO under accession code GSE167665. All data analysed during this study are included in the manuscript and supporting files. Source data files have been provided for all figures and figure supplements. Four additional and previously published datasets are used in this study: Jonker et al. 2013, GTEx Consortium et al. 2017, Schaum et al. 2020, and Tabula Muris Consortium 2020. All the code used to perform analyses is available in GitHub: https://github.com/hmtzg/geneexp_mouse (copy archived at swh:1:rev:1f2434f90404a79c87d545eca8723d99b123ac1c).

1. 1. Izgi H
   2. Han D
   3. Isildak U
   4. Huang S
   5. Kocabiyik E
   6. Khaitovich P
   7. Somel M
   8. Donertas HM

   (2021) NCBI Gene Expression Omnibus

   ID GSE167665. Bulk RNA-seq of mice covering the whole lifespan (2 days to 904 days) from four tissues.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE167665

1. 1. Jonker MJ
   2. Melis JP
   3. Kuiper RV
   4. van der Hoeven TV
   5. Robinson J
   6. van der Horst GT
   7. Breit TM
   8. Vijg J
   9. Dollé ME
   10. Hoeijmakers JH
   11. van Steeg H

   (2013) NCBI Gene Expression Omnibus

   ID GSE34378. Aging Experiment.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34378
2. 1. GTEx Consortium et al

   (2017) GTEx Portal

   ID phs000424.v8.p2. Gene TPMs.

   https://storage.googleapis.com/gtex_analysis_v8/rna_seq_data/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_tpm.gct.gz
3. 1. Pisco A

   (2020) figshare

   Official data release for Tabula Muris Senis.

   https://doi.org/10.6084/m9.figshare.12654728.v1
4. 1. Schaum et al

   (2020) NCBI Gene Expression Omnibus

   ID GSE132040. Tabula Muris Senis: Bulk sequencing.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132040

1. Software

   1. Alexa A
   2. Rahnenfuhrer J

   (2019)

   topGO: Enrichment Analysis for Gene Ontology

   TopGO.
2. 1. Anders S
   2. Huber W

   (2010) Differential expression analysis for sequence count data

   Genome Biology 11:R106.

   https://doi.org/10.1186/gb-2010-11-10-r106

   • PubMed
   • Google Scholar
3. 1. Anders S
   2. Pyl PT
   3. Huber W

   (2015) HTSeq--a Python framework to work with high-throughput sequencing data

   Bioinformatics (Oxford, England) 31:166–169.

   https://doi.org/10.1093/bioinformatics/btu638

   • PubMed
   • Google Scholar
4. Software

   1. Andrews S

   (2010) FastQC: A Quality Control Tool for High Throughput Sequence Data

   FastQC.

   http://www.bioinformatics.babraham.ac.uk/projects/fastqc
5. 1. Angelidis I
   2. Simon LM
   3. Fernandez IE
   4. Strunz M
   5. Mayr CH
   6. Greiffo FR
   7. Tsitsiridis G
   8. Ansari M
   9. Graf E
   10. Strom T-M
   11. Nagendran M
   12. Desai T
   13. Eickelberg O
   14. Mann M
   15. Theis FJ
   16. Schiller HB

   (2019) An atlas of the aging lung mapped by single cell transcriptomics and deep tissue proteomics

   Nature Communications 10:963.

   https://doi.org/10.1038/s41467-019-08831-9

   • PubMed
   • Google Scholar
6. 1. Anisimova AS
   2. Meerson MB
   3. Gerashchenko MV
   4. Kulakovskiy IV
   5. Dmitriev SE
   6. Gladyshev VN

   (2020) Multifaceted deregulation of gene expression and protein synthesis with age

   PNAS 117:15581–15590.

   https://doi.org/10.1073/pnas.2001788117

   • PubMed
   • Google Scholar
7. 1. Ashburner M
   2. Ball CA
   3. Blake JA
   4. Botstein D
   5. Butler H
   6. Cherry JM
   7. Davis AP
   8. Dolinski K
   9. Dwight SS
   10. Eppig JT
   11. Harris MA
   12. Hill DP
   13. Issel-Tarver L
   14. Kasarskis A
   15. Lewis S
   16. Matese JC
   17. Richardson JE
   18. Ringwald M
   19. Rubin GM
   20. Sherlock G

   (2000) Gene Ontology: tool for the unification of biology

   Nature Genetics 25:25–29.

   https://doi.org/10.1038/75556

   • Google Scholar
8. 1. Bahar R
   2. Hartmann CH
   3. Rodriguez KA
   4. Denny AD
   5. Busuttil RA
   6. Dollé MET
   7. Calder RB
   8. Chisholm GB
   9. Pollock BH
   10. Klein CA
   11. Vijg J

   (2006) Increased cell-to-cell variation in gene expression in ageing mouse heart

   Nature 441:1011–1014.

   https://doi.org/10.1038/nature04844

   • PubMed
   • Google Scholar
9. 1. Battle A
   2. Brown CD
   3. Engelhardt BE
   4. Montgomery SB
   5. GTEx Consortium
   6. Laboratory, Data Analysis &Coordinating Center (LDACC)
   7. Statistical Methods groups—Analysis Working Group
   8. Enhancing GTEx (eGTEx) groups
   9. NIH Common Fund
   10. NIH/NCI
   11. NIH/NHGRI
   12. NIH/NIMH
   13. NIH/NIDA
   14. Biospecimen Collection Source Site—NDRI
   15. Biospecimen Collection Source Site—RPCI
   16. Biospecimen Core Resource—VARI
   17. Brain Bank Repository—University of Miami Brain Endowment Bank
   18. Leidos Biomedical—Project Management
   19. ELSI Study
   20. Genome Browser Data Integration &Visualization—EBI
   21. Genome Browser Data Integration &Visualization—UCSC Genomics Institute, University of California Santa Cruz
   22. Lead analysts
   23. NIH program management
   24. Biospecimen collection
   25. Pathology
   26. eQTL manuscript working group

   (2017) Genetic effects on gene expression across human tissues

   Nature 550:204–213.

   https://doi.org/10.1038/nature24277

   • PubMed
   • Google Scholar
10. 1. Benjamini Y
    2. Hochberg Y

    (1995) Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing

    Journal of the Royal Statistical Society 57:289–300.

    https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

    • Google Scholar
11. 1. Blagosklonny MV

    (2006) Aging and immortality: quasi-programmed senescence and its pharmacologic inhibition

    Cell Cycle (Georgetown, Tex.) 5:2087–2102.

    https://doi.org/10.4161/cc.5.18.3288

    • PubMed
    • Google Scholar
12. 1. Bolger AM
    2. Lohse M
    3. Usadel B

    (2014) Trimmomatic: a flexible trimmer for Illumina sequence data

    Bioinformatics (Oxford, England) 30:2114–2120.

    https://doi.org/10.1093/bioinformatics/btu170

    • PubMed
    • Google Scholar
13. Software

    1. Bolstad B

    (2020) preprocessCore: A Collection of Pre-Processing Functions

    PreprocessCore.

    https://github.com/bmbolstad/preprocessCore
14. 1. Brawand D
    2. Soumillon M
    3. Necsulea A
    4. Julien P
    5. Csárdi G
    6. Harrigan P
    7. Weier M
    8. Liechti A
    9. Aximu-Petri A
    10. Kircher M
    11. Albert FW
    12. Zeller U
    13. Khaitovich P
    14. Grützner F
    15. Bergmann S
    16. Nielsen R
    17. Pääbo S
    18. Kaessmann H

    (2011) The evolution of gene expression levels in mammalian organs

    Nature 478:343–348.

    https://doi.org/10.1038/nature10532

    • PubMed
    • Google Scholar
15. 1. Cardoso-Moreira M
    2. Halbert J
    3. Valloton D
    4. Velten B
    5. Chen C
    6. Shao Y
    7. Liechti A
    8. Ascenção K
    9. Rummel C
    10. Ovchinnikova S
    11. Mazin PV
    12. Xenarios I
    13. Harshman K
    14. Mort M
    15. Cooper DN
    16. Sandi C
    17. Soares MJ
    18. Ferreira PG
    19. Afonso S
    20. Carneiro M
    21. Turner JMA
    22. VandeBerg JL
    23. Fallahshahroudi A
    24. Jensen P
    25. Behr R
    26. Lisgo S
    27. Lindsay S
    28. Khaitovich P
    29. Huber W
    30. Baker J
    31. Anders S
    32. Zhang YE
    33. Kaessmann H

    (2019) Gene expression across mammalian organ development

    Nature 571:505–509.

    https://doi.org/10.1038/s41586-019-1338-5

    • PubMed
    • Google Scholar
16. 1. Carvalho BS
    2. Irizarry RA

    (2010) A framework for oligonucleotide microarray preprocessing

    Bioinformatics (Oxford, England) 26:2363–2367.

    https://doi.org/10.1093/bioinformatics/btq431

    • PubMed
    • Google Scholar
17. 1. Colantuoni C
    2. Lipska BK
    3. Ye T
    4. Hyde TM
    5. Tao R
    6. Leek JT
    7. Colantuoni EA
    8. Elkahloun AG
    9. Herman MM
    10. Weinberger DR
    11. Kleinman JE

    (2011) Temporal dynamics and genetic control of transcription in the human prefrontal cortex

    Nature 478:519–523.

    https://doi.org/10.1038/nature10524

    • PubMed
    • Google Scholar
18. 1. de Magalhães JP
    2. Church GM

    (2005) Genomes optimize reproduction: aging as a consequence of the developmental program

    Physiology (Bethesda, Md.) 20:252–259.

    https://doi.org/10.1152/physiol.00010.2005

    • PubMed
    • Google Scholar
19. 1. Dobin A
    2. Davis CA
    3. Schlesinger F
    4. Drenkow J
    5. Zaleski C
    6. Jha S
    7. Batut P
    8. Chaisson M
    9. Gingeras TR

    (2013) STAR: ultrafast universal RNA-seq aligner

    Bioinformatics (Oxford, England) 29:15–21.

    https://doi.org/10.1093/bioinformatics/bts635

    • PubMed
    • Google Scholar
20. 1. Dönertaş HM
    2. İzgi H
    3. Kamacıoğlu A
    4. He Z
    5. Khaitovich P
    6. Somel M

    (2017) Gene expression reversal toward pre-adult levels in the aging human brain and age-related loss of cellular identity

    Scientific Reports 7:5894.

    https://doi.org/10.1038/s41598-017-05927-4

    • PubMed
    • Google Scholar
21. 1. Dönertaş HM
    2. Fuentealba Valenzuela M
    3. Partridge L
    4. Thornton JM

    (2018) Gene expression-based drug repurposing to target aging

    Aging Cell 17:e12819.

    https://doi.org/10.1111/acel.12819

    • PubMed
    • Google Scholar
22. 1. Dönertaş HM
    2. Fabian DK
    3. Valenzuela MF
    4. Partridge L
    5. Thornton JM

    (2021) Common genetic associations between age-related diseases

    Nature Aging 1:400–412.

    https://doi.org/10.1038/s43587-021-00051-5

    • PubMed
    • Google Scholar
23. 1. Durinck S
    2. Spellman PT
    3. Birney E
    4. Huber W

    (2009) Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt

    Nature Protocols 4:1184–1191.

    https://doi.org/10.1038/nprot.2009.97

    • PubMed
    • Google Scholar
24. 1. Enge M
    2. Arda HE
    3. Mignardi M
    4. Beausang J
    5. Bottino R
    6. Kim SK
    7. Quake SR

    (2017) Single-Cell Analysis of Human Pancreas Reveals Transcriptional Signatures of Aging and Somatic Mutation Patterns

    Cell 171:321–330.

    https://doi.org/10.1016/j.cell.2017.09.004

    • PubMed
    • Google Scholar
25. 1. Ezcurra M
    2. Benedetto A
    3. Sornda T
    4. Gilliat AF
    5. Au C
    6. Zhang Q
    7. van Schelt S
    8. Petrache AL
    9. Wang H
    10. Guardia Y
    11. Bar-Nun S
    12. Tyler E
    13. Wakelam MJ
    14. Gems D

    (2018) C. elegans Eats Its Own Intestine to Make Yolk Leading to Multiple Senescent Pathologies

    Current Biology 28:3352.

    https://doi.org/10.1016/j.cub.2018.10.003

    • PubMed
    • Google Scholar
26. 1. Feser J
    2. Truong D
    3. Das C
    4. Carson JJ
    5. Kieft J
    6. Harkness T
    7. Tyler JK

    (2010) Elevated histone expression promotes life span extension

    Molecular Cell 39:724–735.

    https://doi.org/10.1016/j.molcel.2010.08.015

    • PubMed
    • Google Scholar
27. Book

    1. Fisher RA

    (1930) The Genetical Theory of Natural Selection

    OUP.

    https://doi.org/10.5962/bhl.title.27468

    • Google Scholar
28. Book

    1. Flurkey K
    2. Currer JM
    3. Harrison DE

    (2007)

    Chapter 20 - Mouse Models in Aging Research

    In: Fox James G, Davisson Muriel T, Quimby Fred W, Barthold Stephen W, Newcomer Christian E, Smith Abigail L, editors. The Mouse in Biomedical Research. Academic Press. pp. 637–672.

    • Google Scholar
29. Book

    1. Fox J
    2. Weisberg S

    (2018)

    An R Companion to Applied Regression

    SAGE Publications.

    • Google Scholar
30. 1. Gems D
    2. Partridge L

    (2013) Genetics of longevity in model organisms: debates and paradigm shifts

    Annual Review of Physiology 75:621–644.

    https://doi.org/10.1146/annurev-physiol-030212-183712

    • PubMed
    • Google Scholar
31. 1. Hernando-Herraez I
    2. Evano B
    3. Stubbs T
    4. Commere P-H
    5. Jan Bonder M
    6. Clark S
    7. Andrews S
    8. Tajbakhsh S
    9. Reik W

    (2019) Ageing affects DNA methylation drift and transcriptional cell-to-cell variability in mouse muscle stem cells

    Nature Communications 10:4361.

    https://doi.org/10.1038/s41467-019-12293-4

    • Google Scholar
32. Software

    1. hmtzg

    (2022) geneexp_mouse, version swh:1:rev:1f2434f90404a79c87d545eca8723d99b123ac1c

    Software Heritage.

    https://archive.softwareheritage.org/swh:1:dir:b8c43e421f7216167380682c06ed9040db053627;origin=https://github.com/hmtzg/geneexp_mouse;visit=swh:1:snp:5a896cb4722794c85f464a75d459caf84021ffa0;anchor=swh:1:rev:1f2434f90404a79c87d545eca8723d99b123ac1c
33. 1. Hsu SD
    2. Lin FM
    3. Wu WY
    4. Liang C
    5. Huang WC
    6. Chan WL
    7. Tsai WT
    8. Chen GZ
    9. Lee CJ
    10. Chiu CM
    11. Chien CH
    12. Wu MC
    13. Huang CY
    14. Tsou AP
    15. Huang HD

    (2011) miRTarBase: a database curates experimentally validated microRNA-target interactions

    Nucleic Acids Research 39:D163–D169.

    https://doi.org/10.1093/nar/gkq1107

    • PubMed
    • Google Scholar
34. 1. Hsu SD
    2. Tseng YT
    3. Shrestha S
    4. Lin YL
    5. Khaleel A
    6. Chou CH
    7. Chu CF
    8. Huang HY
    9. Lin CM
    10. Ho SY
    11. Jian TY
    12. Lin FM
    13. Chang TH
    14. Weng SL
    15. Liao KW
    16. Liao IE
    17. Liu CC
    18. Huang HD

    (2014) miRTarBase update 2014: an information resource for experimentally validated miRNA-target interactions

    Nucleic Acids Research 42:D78–D85.

    https://doi.org/10.1093/nar/gkt1266

    • PubMed
    • Google Scholar
35. 1. Işıldak U
    2. Somel M
    3. Thornton JM
    4. Dönertaş HM

    (2020) Temporal changes in the gene expression heterogeneity during brain development and aging

    Scientific Reports 10:4080.

    https://doi.org/10.1038/s41598-020-60998-0

    • PubMed
    • Google Scholar
36. 1. Jonker MJ
    2. Melis JP
    3. Kuiper RV
    4. van der Hoeven TV
    5. Wackers PFK
    6. Robinson J
    7. van der Horst GT
    8. Dollé ME
    9. Vijg J
    10. Breit TM
    11. Hoeijmakers JH
    12. van Steeg H

    (2013) Life spanning murine gene expression profiles in relation to chronological and pathological aging in multiple organs

    Aging Cell 12:901–909.

    https://doi.org/10.1111/acel.12118

    • PubMed
    • Google Scholar
37. 1. Kedlian VR
    2. Donertas HM
    3. Thornton JM

    (2019) The widespread increase in inter-individual variability of gene expression in the human brain with age

    Aging 11:2253–2280.

    https://doi.org/10.18632/aging.101912

    • PubMed
    • Google Scholar
38. 1. Kim S
    2. Villeponteau B
    3. Jazwinski SM

    (1996) Effect of replicative age on transcriptional silencing near telomeres in Saccharomyces cerevisiae

    Biochemical and Biophysical Research Communications 219:370–376.

    https://doi.org/10.1006/bbrc.1996.0240

    • PubMed
    • Google Scholar
39. 1. Lind MI
    2. Ravindran S
    3. Sekajova Z
    4. Carlsson H
    5. Hinas A
    6. Maklakov AA

    (2019) Experimentally reduced insulin/IGF-1 signaling in adulthood extends lifespan of parents and improves Darwinian fitness of their offspring

    Evolution Letters 3:207–216.

    https://doi.org/10.1002/evl3.108

    • PubMed
    • Google Scholar
40. 1. Liu X
    2. Han D
    3. Somel M
    4. Jiang X
    5. Hu H
    6. Guijarro P
    7. Zhang N
    8. Mitchell A
    9. Halene T
    10. Ely JJ
    11. Sherwood CC
    12. Hof PR
    13. Qiu Z
    14. Pääbo S
    15. Akbarian S
    16. Khaitovich P

    (2016) Disruption of an Evolutionarily Novel Synaptic Expression Pattern in Autism

    PLOS Biology 14:e1002558.

    https://doi.org/10.1371/journal.pbio.1002558

    • PubMed
    • Google Scholar
41. 1. Love MI
    2. Huber W
    3. Anders S

    (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

    Genome Biology 15:550.

    https://doi.org/10.1186/s13059-014-0550-8

    • PubMed
    • Google Scholar
42. 1. Luegmayr E
    2. Glantschnig H
    3. Wesolowski GA
    4. Gentile MA
    5. Fisher JE
    6. Rodan GA
    7. Reszka AA

    (2004) Osteoclast formation, survival and morphology are highly dependent on exogenous cholesterol/lipoproteins

    Cell Death and Differentiation 11 Suppl 1:S108–S118.

    https://doi.org/10.1038/sj.cdd.4401399

    • PubMed
    • Google Scholar
43. 1. Martinez-Jimenez CP
    2. Eling N
    3. Chen H-C
    4. Vallejos CA
    5. Kolodziejczyk AA
    6. Connor F
    7. Stojic L
    8. Rayner TF
    9. Stubbington MJT
    10. Teichmann SA
    11. de la Roche M
    12. Marioni JC
    13. Odom DT

    (2017) Aging increases cell-to-cell transcriptional variability upon immune stimulation

    Science (New York, N.Y.) 355:1433–1436.

    https://doi.org/10.1126/science.aah4115

    • PubMed
    • Google Scholar
44. 1. Matys V
    2. Fricke E
    3. Geffers R
    4. Gössling E
    5. Haubrock M
    6. Hehl R
    7. Hornischer K
    8. Karas D
    9. Kel AE
    10. Kel-Margoulis OV
    11. Kloos D-U
    12. Land S
    13. Lewicki-Potapov B
    14. Michael H
    15. Münch R
    16. Reuter I
    17. Rotert S
    18. Saxel H
    19. Scheer M
    20. Thiele S
    21. Wingender E

    (2003) TRANSFAC: transcriptional regulation, from patterns to profiles

    Nucleic Acids Research 31:374–378.

    https://doi.org/10.1093/nar/gkg108

    • PubMed
    • Google Scholar
45. 1. Matys V
    2. Kel-Margoulis OV
    3. Fricke E
    4. Liebich I
    5. Land S
    6. Barre-Dirrie A
    7. Reuter I
    8. Chekmenev D
    9. Krull M
    10. Hornischer K
    11. Voss N
    12. Stegmaier P
    13. Lewicki-Potapov B
    14. Saxel H
    15. Kel AE
    16. Wingender E

    (2006) TRANSFAC and its module TRANSCompel: transcriptional gene regulation in eukaryotes

    Nucleic Acids Research 34:D108–D110.

    https://doi.org/10.1093/nar/gkj143

    • PubMed
    • Google Scholar
46. 1. Medawar PB

    (1953)

    Unsolved problem of biology

    The Medical Journal of Australia 1:854–855.

    • Google Scholar
47. Software

    1. Pisco A

    (2020) Tabula Muris Senis Data Objects

    Figshare.

    https://doi.org/10.6084/m9.figshare.12654728.v1
48. Software

    1. Raivo K

    (2019)

    Pheatmap: Pretty Heatmaps

    R Package Version.
49. 1. Rouillard AD
    2. Gundersen GW
    3. Fernandez NF
    4. Wang Z
    5. Monteiro CD
    6. McDermott MG
    7. Ma’ayan A

    (2016) The harmonizome: a collection of processed datasets gathered to serve and mine knowledge about genes and proteins

    Database 2016:baw100.

    https://doi.org/10.1093/database/baw100

    • PubMed
    • Google Scholar
50. 1. Sampathkumar NK
    2. Bravo JI
    3. Chen Y
    4. Danthi PS
    5. Donahue EK
    6. Lai RW
    7. Lu R
    8. Randall LT
    9. Vinson N
    10. Benayoun BA

    (2020) Widespread sex dimorphism in aging and age-related diseases

    Human Genetics 139:333–356.

    https://doi.org/10.1007/s00439-019-02082-w

    • PubMed
    • Google Scholar
51. 1. Schaum N
    2. Lehallier B
    3. Hahn O
    4. Pálovics R
    5. Hosseinzadeh S
    6. Lee SE
    7. Sit R
    8. Lee DP
    9. Losada PM
    10. Zardeneta ME
    11. Fehlmann T
    12. Webber JT
    13. McGeever A
    14. Calcuttawala K
    15. Zhang H
    16. Berdnik D
    17. Mathur V
    18. Tan W
    19. Zee A
    20. Tan M
    21. Tabula Muris Consortium
    22. Pisco AO
    23. Karkanias J
    24. Neff NF
    25. Keller A
    26. Darmanis S
    27. Quake SR
    28. Wyss-Coray T

    (2020) Ageing hallmarks exhibit organ-specific temporal signatures

    Nature 583:596–602.

    https://doi.org/10.1038/s41586-020-2499-y

    • PubMed
    • Google Scholar
52. 1. Somel M
    2. Khaitovich P
    3. Bahn S
    4. Pääbo S
    5. Lachmann M

    (2006) Gene expression becomes heterogeneous with age

    Current Biology 16:R359–R360.

    https://doi.org/10.1016/j.cub.2006.04.024

    • PubMed
    • Google Scholar
53. 1. Somel M
    2. Guo S
    3. Fu N
    4. Yan Z
    5. Hu HY
    6. Xu Y
    7. Yuan Y
    8. Ning Z
    9. Hu Y
    10. Menzel C
    11. Hu H
    12. Lachmann M
    13. Zeng R
    14. Chen W
    15. Khaitovich P

    (2010) MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain

    Genome Research 20:1207–1218.

    https://doi.org/10.1101/gr.106849.110

    • PubMed
    • Google Scholar
54. 1. Stuart T
    2. Butler A
    3. Hoffman P
    4. Hafemeister C
    5. Papalexi E
    6. Mauck WM
    7. Hao Y
    8. Stoeckius M
    9. Smibert P
    10. Satija R

    (2019) Comprehensive Integration of Single-Cell Data

    Cell 177:1888–1902.

    https://doi.org/10.1016/j.cell.2019.05.031

    • PubMed
    • Google Scholar
55. 1. Tabula Muris Consortium

    (2020) A single-cell transcriptomic atlas characterizes ageing tissues in the mouse

    Nature 583:590–595.

    https://doi.org/10.1038/s41586-020-2496-1

    • PubMed
    • Google Scholar
56. 1. Tacutu R
    2. Thornton D
    3. Johnson E
    4. Budovsky A
    5. Barardo D
    6. Craig T
    7. Diana E
    8. Lehmann G
    9. Toren D
    10. Wang J
    11. Fraifeld VE
    12. de Magalhães JP

    (2018) Human Ageing Genomic Resources: new and updated databases

    Nucleic Acids Research 46:D1083–D1090.

    https://doi.org/10.1093/nar/gkx1042

    • PubMed
    • Google Scholar
57. 1. Tibshirani R
    2. Walther G
    3. Hastie T

    (2001) Estimating the number of clusters in a data set via the gap statistic

    Journal of the Royal Statistical Society 63:411–423.

    https://doi.org/10.1111/1467-9868.00293

    • Google Scholar
58. 1. Trapnell C
    2. Williams BA
    3. Pertea G
    4. Mortazavi A
    5. Kwan G
    6. van Baren MJ
    7. Salzberg SL
    8. Wold BJ
    9. Pachter L

    (2010) Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation

    Nature Biotechnology 28:511–515.

    https://doi.org/10.1038/nbt.1621

    • PubMed
    • Google Scholar
59. 1. Turan ZG
    2. Parvizi P
    3. Dönertaş HM
    4. Tung J
    5. Khaitovich P
    6. Somel M

    (2019) Molecular footprint of Medawar’s mutation accumulation process in mammalian aging

    Aging Cell 18:e12965.

    https://doi.org/10.1111/acel.12965

    • PubMed
    • Google Scholar
60. 1. Williams GC

    (1957) PLEIOTROPY, NATURAL SELECTION, AND THE EVOLUTION OF SENESCENCE

    Evolution 11:398–411.

    https://doi.org/10.1111/j.1558-5646.1957.tb02911.x

    • Google Scholar
61. 1. Yang J
    2. Huang T
    3. Petralia F
    4. Long Q
    5. Zhang B
    6. Argmann C
    7. Zhao Y
    8. Mobbs CV
    9. Schadt EE
    10. Zhu J
    11. Tu Z
    12. GTEx Consortium

    (2015) Synchronized age-related gene expression changes across multiple tissues in human and the link to complex diseases

    Scientific Reports 5:15145.

    https://doi.org/10.1038/srep15145

    • PubMed
    • Google Scholar
62. Book

    1. Yang J
    2. Griffin P
    3. Vera D
    4. Apostolides J
    5. Hayano M
    6. Meer M
    7. Salfati E

    (2019)

    Erosion of the Epigenetic Landscape and Loss of Cellular Identity as a Cause of Aging in Mammals

    Cold Spring Harbor Laboratory.

    • Google Scholar
63. 1. Yu G
    2. Wang L-G
    3. Han Y
    4. He Q-Y

    (2012) clusterProfiler: an R package for comparing biological themes among gene clusters

    Omics 16:284–287.

    https://doi.org/10.1089/omi.2011.0118

    • PubMed
    • Google Scholar
64. 1. Yuan Y
    2. Chen Y-PP
    3. Boyd-Kirkup J
    4. Khaitovich P
    5. Somel M

    (2012) Accelerated aging-related transcriptome changes in the female prefrontal cortex

    Aging Cell 11:894–901.

    https://doi.org/10.1111/j.1474-9726.2012.00859.x

    • PubMed
    • Google Scholar
65. 1. Zahn JM
    2. Poosala S
    3. Owen AB
    4. Ingram DK
    5. Lustig A
    6. Carter A
    7. Weeraratna AT
    8. Taub DD
    9. Gorospe M
    10. Mazan-Mamczarz K
    11. Lakatta EG
    12. Boheler KR
    13. Xu X
    14. Mattson MP
    15. Falco G
    16. Ko MSH
    17. Schlessinger D
    18. Firman J
    19. Kummerfeld SK
    20. Wood WH
    21. Zonderman AB
    22. Kim SK
    23. Becker KG

    (2007) AGEMAP: a gene expression database for aging in mice

    PLOS Genetics 3:e201.

    https://doi.org/10.1371/journal.pgen.0030201

    • PubMed
    • Google Scholar

Decision letter after peer review:

Thank you for submitting your article "Inter-tissue convergence of gene expression and loss of cellular identity during ageing" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Matt Kaeberlein as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this letter to help you prepare a revised submission.

Essential revisions:

1. The reviewers noted some potential issues with the normalization of the RNA-seq data that could affect the results, including the use of FPKM normalization, which may drive some of the observations due to lack of inter-sample-normalization (see specific comments by Reviewers #2 and 3). The authors would need to show that using appropriate count-based methods yields similar conclusions.

2. The authors often support their conclusions on low effect sizes and significance above usual thresholds. It would be important to explicitly discuss these as caveats, as well as strongly tone down these conclusions (including in the title.). More generally, the reviewers felt that potential alternative hypotheses needed to be clearly laid out, and potential ruled out, for the manuscript to stand more strongly (see specific comments by Reviewers #1, 2 and 3.)

3. In general, the reviewers felt that additional information on methods was required for reproducibility (see individual reviewer reports for specific points).

Reviewer #1 (Recommendations for the authors):

1. "Such "runaway development" phenomena may be due to insufficient natural selection pressure to optimise regulation after sexual reproduction starts." As stated, this sentence suggests that maladaptive traits that could lower an organism's total reproductive success would not be selected against as long as they only appear after the beginning of reproduction. I would perhaps ask that here the authors cite, not a general review of evolutionary theories of aging, but a publication that specifically proposes that alleles that are deleterious as early as the beginning of sexual reproduction nevertheless are not selected against.

2. Lines 84-91, 145-153, 229-240, 245-250, etc. Here the font size changes in the manuscript which I assume is not intended.

3. Figure 1 caption. "Similarities were calculated using Spearman correlations coefficient between expression-age correlations across tissues." Here I suspect that "Spearman correlations coefficient" was intended to be "Spearman's correlation coefficient".

Reviewer #2 (Recommendations for the authors):

Although the authors report an intriguing finding, there are major issues in the manuscript as it stands, notably concerning the clarity and rigor of the data analysis and manuscript.

1. For the analysis of their RNA-seq dataset across samples, the authors describe using the FPKM framework (methods, line 438). This is very problematic for a number of reasons – use of FPKM for differential gene analysis across samples is deprecated because it poorly controls for inter-sample variability [PMID: 22988256; PMID: 32284352]. FPKM are not comparable between samples because the normalization is sample-specific – thus, since no cross-library normalization step is performed, it should never be used for between sample comparison. Thanks to community benchmarking efforts, TMM or VST normalized counts are widely believed to be the superior choice for DGE analysis [PMID: 22988256; PMID: 32284352]. Thus, it is crucial that the authors correct their analysis to avoid the use of FPKM, and show that their conclusions would hold with an analysis on TMM or VST normalized counts.

2. There are issues regarding consistency in analytical choices throughout the paper.

a. Application of significance cutoffs vary widely in the manuscript when selecting specific gene sets. Reasoning should be provided for each case in which significance cutoff was (or not) applied (e.g. Line 111). For example, it is rather concerning that in Figure 1d, no gene was significantly up-regulated with age in the muscle (significance cutoff applied), but in Figure 1f, a large proportion of genes are shown to increase in expression with age in the muscle (significance cutoff not applied).

b. Authors are not consistent with the FDR threshold applied (e.g. FDR < 10% in Line 123, FDR 20% in Line 127). The reasoning for the chosen thresholds, as well as the different choices should be provided for cases in which higher than standard FDR is applied.

c. In Figure 1, authors show that the muscle presents with the least number of significantly changed gene expression patterns (up and down) with age. Additionally, in the manuscript, the authors explain that the lack of genes that show statistically significant change in the muscle is supported by a previous study that also showed "weak ageing transcriptome signature across multiple datasets (Line 137)." Thus, it is rather concerning that excluding the cortex, a tissue with a large number of differentially expressed genes, for CoV calculation resulted in greater significance.

d. The authors use a mix of multiple correction methods without clear explanations as to why a consistent measure isn't used (i.e. FDR/BH, BY) – BY test is justified for GO analysis line 472, but why BY was not used for DGE analysis is not clear (line 462).

e. A number of statistical analyses performed in the manuscript do not pass the standard significance thresholds (FDR < 5% or p < 0.05). The standard applied in the field is usually p-value or FDR < 5%, although relaxing this threshold can be done when dealing with noisier datasets. When using more relaxed FDR or p-value thresholds, the authors should also be careful to strongly tone down statements proportionally to the lower statistical support for their claims (eg. Line 240: "a SIGNIFICANT divergence-convergence pattern).

3. The authors use a dataset from exclusively male mice. Since aging is known to be extremely sex-dimorphic [PMID: 27304504], the authors need to explicitly discuss this caveat/limitation in the main text (i.e. these observations may not hold in female animals). If the authors wish to make a more general statement, a suggestion would be to reanalyze data from the recent Tabula Muris Senis Companion paper Schaum et al., [PMID: 32669715], which includes similar tissues, "developmental" time points (1 vs 3 months) and aging course (3 to 27 months), and both sexes to extend the findings. In addition, the Schaum dataset is also RNA-seq, and thus more directly comparable to the author's dataset compared to the microarray Jonker dataset used for partial validation (Figure 2 – supplement 8).

4. The authors should consider revising the title, since the current form seems to indicate that the work discovered s significant loss of cellular identity during aging. This statement is too strong based on the analyses of the paper. We suggest toning down the title, for instance to "Transcriptomic analyses suggest inter-tissue convergence of gene expression and loss of cellular identity during ageing".

Reviewer #3 (Recommendations for the authors):

Authors only analyze the gene sets that are correlated with age i.e. linearly go up or down with age. However, ageing trajectories for different genes can be very dynamic. This is the likely reason why only a few age-related changes were detected in some tissues. To assess the robustness of some of the key observations, authors should consider an orthogonal and unbiased clustering analysis similar to Figure 2 supplement to get the patterns for gene expression changes with age and to what extent DiCo is observed with different clusters.

Some of the key observations have either low effect size or are not statistically significant. For example, inter-tissue expression similarity during ageing and development (lines 115-119), the expression reversal (lines 158-161), changes in coefficient of variation with age (lines 202-207) etc. Authors clearly mention these in the manuscript, but it raises concerns on the extent to which some of these observations might be random, or due to idiosyncrasies of RNA-seq itself. These caveats should also be discussed.

In the global PCA based analysis for figure 1 and 2 combining both development and ageing data, the variance can be driven by one of the two states. For example, transcription divergence analysis in PCA in figure 1 based on development alone is not significant which suggests that the divergence of the developmental samples is driven in-fact by ageing. For PCA based analysis using separate PCAs for development and ageing will be more statistically sound.

Are there any known gene categories that have stronger expression reversal trends? For example, do known transcription factors or pathways that regulate development show DiCo if analyzed separately? Or is this trend only seen at the genome-wide level and completely stochastic?

For functional enrichment analyses authors only mention some selected functions in the text but GO terms have a lot of redundancy and it is unclear if there is any emerging functional trends looking beyond the redundancy or selected examples. It will be useful to thoroughly analyze and summarize the enriched GO terms to assess for any functional patterns that connect functional decline with ageing and convergent expression or DiCo.

Is there a significant overlap between mouse and human ageing related changes at the gene or pathway level?

Line 299 – linking the functional categories with cellular identity appears a bit hand-wavy. Stronger support would be needed to make this conclusion.

Some of the aspects need more methodological details:

For RNA-seq analysis, were all the samples normalized together, or were development and ageing samples normalized separately? It should be clarified in the methods.

Were the mouse perfused? The details should be added in methods. If the mouse were not perfused, a significant shared proportion of gene expression can be due to blood components. Also, the possibility that some of the shared changes could be due to infiltration of immune and other cell types in tissues should be mentioned.

More details for single cell deconvolution should be added in methods.

It is unclear to me how how the age related gene expression and PC correlation was performed. Was the loadings underlying PC axes used for correlation analysis? It should be elaborated in methods.

Some of the supplemental figures are missing.

Essential revisions:

1. The reviewers noted some potential issues with the normalization of the RNA-seq data that could affect the results, including the use of FPKM normalization, which may drive some of the observations due to lack of inter-sample-normalization (see specific comments by Reviewers #2 and 3). The authors would need to show that using appropriate count-based methods yields similar conclusions.

We believe this comment mainly stems from the lack of clarity of our methods, which we improved in this version. We had used quantile normalisation on FPKM normalised values, which is an inter-sample-normalisation and normalises for overall library sizes between samples. In this revised version, we further repeat the analysis with VST normalisation following the DESeq2 pipeline and confirm the robustness of our results. The details are provided under the relevant comments below.

2. The authors often support their conclusions on low effect sizes and significance above usual thresholds. It would be important to explicitly discuss these as caveats, as well as strongly tone down these conclusions (including in the title.). More generally, the reviewers felt that potential alternative hypotheses needed to be clearly laid out, and potential ruled out, for the manuscript to stand more strongly (see specific comments by Reviewers #1, 2 and 3.)

We thank the reviewers for raising this issue. We have now (1) explained the difference between FPR (now called eFPP to prevent confusion) and FDR in the methods, to clarify that we use the same significance cutoff across different analyses, and also provide an additional measure that estimates the false positive proportion (not a p-value) based on permutations; (2) toned down the manuscript, including the title; (3) included additional discussion to summarise all the limitations laid out throughout the manuscript; (4) specifically tested the alternative hypothesis suggested by the reviewer #1; (5) proposed ageing-related increase in inter-cellular noise (heteroskedasticity) as a possible explanation for the convergence observed during ageing.

3. In general, the reviewers felt that additional information on methods was required for reproducibility (see individual reviewer reports for specific points).

We have re-written a significant proportion of the methods and included detailed explanations in the figure legends where applicable to improve reproducibility.

Reviewer #1 (Recommendations for the authors):

1. "Such "runaway development" phenomena may be due to insufficient natural selection pressure to optimise regulation after sexual reproduction starts." As stated, this sentence suggests that maladaptive traits that could lower an organism's total reproductive success would not be selected against as long as they only appear after the beginning of reproduction. I would perhaps ask that here the authors cite, not a general review of evolutionary theories of aging, but a publication that specifically proposes that alleles that are deleterious as early as the beginning of sexual reproduction nevertheless are not selected against.

We thank the reviewer for the suggestion. The article by de Magalhaes and Church is the review article in which the authors revised the developmental theory of ageing especially in the light of genomic studies and thus is the original article where the ‘runaway development’ phenomenon was introduced in these exact terms. We now also included the original article by Mikhail V. Blagosklonny, where a related phenomenon, the hyperfunction theory was first proposed (Blagosklonny, 2006) – however, this is again a conceptual review where the theory was first proposed based on observations on many different grounds. Therefore, we have now also included papers with experimental support for this notion (Ezcurra et al., 2018; Lind et al., 2019). Moreover, we have now included a statement in the discussion to convey the message that we observe a reversal in the majority of the genes, but this does not necessarily negate the developmental or hyperfunction theory, which might be relevant for a limited number of pathways and may contribute to the ageing phenotype (lines 550-553). Finally, we now cite (Fisher, 1930; Medawar, 1952; Williams, 1957), which suggests the force of natural selection declines with age.

2. Lines 84-91, 145-153, 229-240, 245-250, etc. Here the font size changes in the manuscript which I assume is not intended.

Thank you for pointing this out – we have now corrected it.

3. Figure 1 caption. "Similarities were calculated using Spearman correlations coefficient between expression-age correlations across tissues." Here I suspect that "Spearman correlations coefficient" was intended to be "Spearman's correlation coefficient".

Thank you – we have changed it and carefully checked the manuscript one more time to prevent such typos.

Reviewer #2 (Recommendations for the authors):

Although the authors report an intriguing finding, there are major issues in the manuscript as it stands, notably concerning the clarity and rigor of the data analysis and manuscript.

1. For the analysis of their RNA-seq dataset across samples, the authors describe using the FPKM framework (methods, line 438). This is very problematic for a number of reasons – use of FPKM for differential gene analysis across samples is deprecated because it poorly controls for inter-sample variability [PMID: 22988256; PMID: 32284352]. FPKM are not comparable between samples because the normalization is sample-specific – thus, since no cross-library normalization step is performed, it should never be used for between sample comparison. Thanks to community benchmarking efforts, TMM or VST normalized counts are widely believed to be the superior choice for DGE analysis [PMID: 22988256; PMID: 32284352]. Thus, it is crucial that the authors correct their analysis to avoid the use of FPKM, and show that their conclusions would hold with an analysis on TMM or VST normalized counts.

We thank the reviewer for their comment. We had used FPKM followed by quantile normalisation, which performs an inter-sample normalisation. We had mentioned this at lines 443, 662, and 682 in the previous version, but had not explained this normalisation procedure in detail. Now the relevant section in the methods is updated to explain the quantile normalisation method (line 690). The main reason we choose quantile normalisation is that we analysed both RNA-seq and microarray data and this normalisation method can be applied to both, following other background normalisations and transformations (i.e. FPKM + log2 transformation for RNA-seq and RMA correction and log2 transformation for microarray). However, we agree with the reviewer that use of VST and TMM are more conventional as these are already implemented in the widely used analysis pipelines of DESeq2 and EdgeR.

We thus repeated the analysis of our dataset using VST and confirmed the main results.

We replicated our findings presented in Figure 1 using VST-normalised data. The results are presented in Figure 1—figure supplement 10 and 11:

i) We first compared VST- and QN-processed data in terms of age-related expression change trends in each tissue and period separately, without applying a significance cutoff (across n=14,973 genes). We found high correlations in all 8 comparisons (⍴_dev=[0.81, 0.94], ⍴_ageing=[0.83, 0.93], Figure 1—figure supplement 11).

ii) Repeating PCA analysis with VST-normalised data, we show that variation in gene expression (n=14,973 genes) is largely explained by tissue differences across PC1, 2, and 3 (ANOVA p<10^-16 for PC1-3). In PC4, which explains 7% of the total variation (as opposed to 8% in Figure 1b, we observe an age-effect across tissues in development (Spearman’s correlation coefficient ⍴=[0.57, 0.99], nominal p<0.01 for each test except in muscle; Figure 1—figure supplement 10b). We also confirmed that tissues diverge during postnatal development, although the convergence during ageing was not significant (change in mean Euclidean distance among tissues with age in PC1-4 space, ⍴_dev=0.937, p_dev=0.00185; during ageing ⍴_ageing=-0.58, p_ageing=0.102, Figure 1-source data).

iii) Similar to Figure 1c results, we find that tissue pairs show weak correlations in their age-related expression change trends measured with VST-normalised data (n=14,973 genes), both during development (⍴=[0.19, 0.39]) and during ageing (⍴=[0.25, 0.34]). Furthermore, we find that the number of genes with the same direction of change (with or without significance cutoff) across four tissues with VST-normalised data is comparable to Figure 1d and Figure 1e result (Figure 1—figure supplement 10d-e).

iv) Reversal analysis with VST-normalised data revealed that ~50% (42-59%) of expressed genes (n=14,973) showed reversal patterns (without significance cutoff) in tissues comparable to Figure 1f result (Figure 1—figure supplement 10f). Comparing VST- and QN-normalised data, we found that 70-79% of the genes showed similar reversal or continuous patterns across tissues.

We confirmed our findings from Figure 2 with VST normalisation and presented the results in Figure 2—figure supplement 14:

i) We calculated CoV values across all expressed genes (n=14,973) among tissues in VST-normalised data and confirmed a significant mean CoV increase in development (Spearman’s correlation coefficient ⍴=0.99, p=10^-5), consistent with QN normalised data that yielded ⍴=0.77 (p=0.041). During ageing, the decrease in mean CoV with age (⍴=-0.48, p=0.23) was comparable to the results shown in Figure 2a (⍴=-0.50, p=0.204).

ii) We replicated the difference between development and ageing in the proportion of genes showing divergent versus convergent patterns (at gene-wise FDR<0.1) with the VSTnormalised data (Figure 2—figure supplement 14d-e). Specifically, VST-normalised data displayed 89% divergence in development (among 3,476 significant genes) as opposed to 70% divergence in QN normalised data (among 2,581 significant genes). Likewise, during ageing, VST-normalised data displayed 68% convergence (among 19 significant genes) similar to 68% convergence in QN normalised data (among 62 significant genes).

We conclude that the choice of the normalisation method does not affect our results. We now present the VST-based results briefly in the text (lines 143-146, 263-265) to prevent redundancy and provide detailed explanations in the supplementary figure legends (Figure 1figure supplement 10, Figure 2—figure supplement 14). We now also explain the VST normalisation steps in the Methods (line 675).

2. There are issues regarding consistency in analytical choices throughout the paper.

a. Application of significance cutoffs vary widely in the manuscript when selecting specific gene sets. Reasoning should be provided for each case in which significance cutoff was (or not) applied (e.g. Line 111). For example, it is rather concerning that in Figure 1d, no gene was significantly up-regulated with age in the muscle (significance cutoff applied), but in Figure 1f, a large proportion of genes are shown to increase in expression with age in the muscle (significance cutoff not applied).

We thank the reviewer for their comment. We actually analyse all the results in two ways: (1) using a consistent significance cutoff (multiple testing (FDR) adjusted p-value<0.1, or if no multiple testing was applied, p-value<0.05), and (2) without using a significance cutoff to assess transcriptome-wide trends. We emphasise this as soon as the first sentence of the Results section titled “Tissues involve common gene expression changes with age”. We now added a second sentence to emphasise this point:

“We next characterised age-related changes in gene expression shared across tissues by (i) studying overall trends at the whole transcriptome level and testing their consistency using permutation tests, and (ii) studying statistically significant changes at the single gene level.”

This has two motivations. First, given the low signal/noise ratios during mammalian ageing, identifying statistically significant patterns in individual genes is frequently difficult, but transcriptome-wide patterns can be more effectively identified. Second, observing the same trends using both approaches increases the robustness of our conclusions and may suggest that the observations are a genome-wide phenomenon, not restricted to a small number of genes. However, since transcriptome-wide trends are more prone to technical artefacts, we further test the genes that show statistically significant changes.

We now attempted to further clarify where a significance cutoff was and was not used in the analyses, by updating the text to emphasise the gene set of interest for each analysis (lines 131, 140).

We also thank the reviewer for pointing out the ambiguity about muscle genes. In Figure 1f, we study global reversal trends in each tissue for all the genes without any significance cutoff, i.e. only analysing expression change trends (following (Dönertaş et al., 2017)). Thus, although there are no statistically significant ageing-related genes in muscle tissue (Figure 1d), we observe a transcriptome-wide reversal trend (Figure 1f). We now updated the figure text to explain clearly what is presented (lines 183-187).

b. Authors are not consistent with the FDR threshold applied (e.g. FDR < 10% in Line 123, FDR 20% in Line 127). The reasoning for the chosen thresholds, as well as the different choices should be provided for cases in which higher than standard FDR is applied.

We thank the reviewer for pointing out that the difference between FDR (e.g. line 123 in the previous version) and FPR (e.g. line 127 in the previous version) was not clearly explained in the text. FDR is the multiple testing correction ‘false discovery rate’ applied using the ‘p.adjust’ function with BH procedure in R, which we did not clearly state in the first submission. We updated the text to explain that FDR is applied to correct the p-values for all statistical tests (functional association tests were corrected with BY procedure in the previous version but for consistency, we used FDR in all the tests in the current version – see point d.) (lines 711-714).

We have also updated the text to replace all “FDR<0.1” with “FDR corrected p-value<0.1”.

We had used the term FPR as an estimate of the false-positive rate calculated using random permutations (FPR = median number of positives observed across permutations / number of positives observed in the data). We provide this metric with permutation test results as additional information on the effect size, not as an alternative to FDR corrected p-values. However, thanks to the reviewer’s comment we realised that the phrase FPR could be confused with other metrics. We, therefore, redefined it as eFPP, ‘estimated false-positive proportion’ (lines 740-743). We now updated the relevant sections to make sure that these two definitions are not interchangeable. We also removed eFPP from the main text to avoid confusion and only provide this information only in the supplementary figures.

c. In Figure 1, authors show that the muscle presents with the least number of significantly changed gene expression patterns (up and down) with age. Additionally, in the manuscript, the authors explain that the lack of genes that show statistically significant change in the muscle is supported by a previous study that also showed "weak ageing transcriptome signature across multiple datasets (Line 137)." Thus, it is rather concerning that excluding the cortex, a tissue with a large number of differentially expressed genes, for CoV calculation resulted in greater significance.

We believe that the reviewer’s concern perhaps stems from our poor wording of CoV calculation in the text. For each individual, we calculated CoV across tissues using all expressed genes (n=15,063) regardless of their specific age-related expression changes (significant or not) in any tissue. Then, we summarised these gene-wise CoV values by calculating their mean to assess a genome-wide variation score among tissues of an individual and how these scores change with age (Figure 2a). We have now updated the relevant Results section to explain the CoV calculation in more detail (lines 221-223). We believe that the significant age-related expression changes in one tissue should not have a direct impact on the results as the number of significantly changed genes during ageing in the cortex was also very small (n=68).

Despite this, we repeated the same analysis, this time excluding the muscle tissue, which resulted in a nonsignificant but similar pattern during ageing (Spearman’s correlation test between mean CoV and age (n=15): ⍴_dev=0.85, p_dev=0.016; ⍴_ageing=-0.29, p_ageing=0.49, Author response image 1a -right panel), and the lack of significance again might be driven by the lack of one cortex sample.

The exclusion of the cortex data was motivated by the fact that we were missing an aged individual in this tissue. However, estimating CoV using only three points may be suboptimal, and we, therefore, decided to be conservative and thus removed this analysis (i.e. results excluding cortex) from the manuscript.

Author response image 1

Download asset Open asset

Change in CoV using 15 samples from three tissues (one oldest individual (904 days of age) excluded due to its missing data in the cortex), either excluding cortex (left column) or muscle (right column) tissue. Each point represents the (a) mean CoV, (b) median CoV value of all protein-coding genes (15,063) for each mouse. x-axis is in log2 scale. The dashed grey line shows the start of the ageing period. The Spearman’s correlation coefficients and p values for each period are indicated separately on the plots.

d. The authors use a mix of multiple correction methods without clear explanations as to why a consistent measure isn't used (i.e. FDR/BH, BY) – BY test is justified for GO analysis line 472, but why BY was not used for DGE analysis is not clear (line 462).

We thank the reviewer for this point. For consistency, we now use FDR (BH) for multiple testing correction for all statistical tests throughout the text. We now updated the relevant sections (lines 711, 713-714).

e. A number of statistical analyses performed in the manuscript do not pass the standard significance thresholds (FDR < 5% or p < 0.05). The standard applied in the field is usually p-value or FDR < 5%, although relaxing this threshold can be done when dealing with noisier datasets. When using more relaxed FDR or p-value thresholds, the authors should also be careful to strongly tone down statements proportionally to the lower statistical support for their claims (eg. Line 240: "a SIGNIFICANT divergence-convergence pattern).

We thank the reviewer for the comment. We believe that the 10% FDR cutoff allows a compromise between overall false positives and false negatives, and is commonly used as a threshold in the transcriptomics literature, including recent analyses on the association between gene expression and splicing (DOI: 10.7554/eLife.67077), expression and neurodegeneration (DOI: 10.7554/eLife.64564), or expression and psychosocial experiences (DOI: 10.7554/eLife.63852).

Using a more stringent threshold would reduce the false positive probability per gene. However, our emphasis in this manuscript is not individual genes but transcriptome-wide patterns, which we test using the whole transcriptome without significance threshold, and also with significant genes (with an FDR corrected p-value<0.1). The particular instance where we use ‘significant’ in line 240 refers to the result of Spearman’s correlation test between mean pairwise correlations among tissues and age in development and ageing periods.

That said, we did follow the reviewer’s suggestion and toned down our conclusions. Accordingly, we changed the referred line to ‘we observed the same divergence-convergence pattern’ (lines 286-287).

3. The authors use a dataset from exclusively male mice. Since aging is known to be extremely sex-dimorphic [PMID: 27304504], the authors need to explicitly discuss this caveat/limitation in the main text (i.e. these observations may not hold in female animals). If the authors wish to make a more general statement, a suggestion would be to reanalyze data from the recent Tabula Muris Senis Companion paper Schaum et al., [PMID: 32669715], which includes similar tissues, "developmental" time points (1 vs 3 months) and aging course (3 to 27 months), and both sexes to extend the findings. In addition, the Schaum dataset is also RNA-seq, and thus more directly comparable to the author's dataset compared to the microarray Jonker dataset used for partial validation (Figure 2 – supplement 8).

We thank the reviewer for this valuable comment. In fact, we identify convergent trends in ageing in both female and male animals, in our dataset (only males) and in the Jonker et al., dataset (only females), respectively. We further find that convergent genes overlap between the two datasets, although modestly (OR = 1.22, p = 1.67x10^-8, Figure 4—figure supplement 2b). This suggests that DiCo is not unique to a single-sex in mice. We now explain this in lines 430-434, 614-616.

We would also have been interested in testing possible sex dimorphism in DiCo, but this is not possible given technical differences between the two datasets (i.e. Jonker et al.,'s and ours).

The GTEx dataset includes both male and female humans and displays patterns consistent with convergence during ageing, although not significantly (Figure 2—figure supplement 811). Following the reviewer's point, we repeated the GTEx dataset analysis separately for females and males. We observed a stronger convergence in females than in males, although none of these results were significant (⍴_female=-0.58, p_female=0.059; ⍴_male= -0.052, p_male=0.77) (Figure 2—figure supplement 16). The female and male GTEx samples differ both in sample size (n=11 vs n=36, respectively) and age range (no male individuals within 20-29 and 70-79 age groups), which could explain the observed differences between sexes.

We now discuss possible sex-dimorphism in Discussion as a limitation of this study and suggest a more comprehensive analysis with different datasets in the future (line 612-620).

Finally, we thank the reviewer for suggesting including the Schaum et al., dataset. We now included this dataset, analysing both the same 4 tissues we have and a broader set of 8 tissues (Figure 2—figure supplement 17-20). However, we could not include the development period as the sample size was very small for the two time points; 1 and 3 months of age. We could not confirm the transcriptome-wide trend towards convergence during ageing in this dataset and found that the muscle and subcutaneous fat tissues showed the most divergent patterns.

However, we confirmed a strong association between the loss of identity and convergence during ageing. Not having the developmental period for the external datasets limited the complete replication of our original study, which focuses on the DiCo pattern, divergence followed by convergence. Nevertheless, a trend towards convergence during ageing was observed at the transcriptome level in Jonker et al., and GTEx and the association between convergence and loss of identity was evident across all datasets. We now updated the text to reflect that the support for the transcriptome-wide trend is weak and our most robust result is the association between identity loss and convergence.

We now add Figure 4—figure supplement 2 and Table 1 to summarise all results across datasets.

4. The authors should consider revising the title, since the current form seems to indicate that the work discovered s significant loss of cellular identity during aging. This statement is too strong based on the analyses of the paper. We suggest toning down the title, for instance to "Transcriptomic analyses suggest inter-tissue convergence of gene expression and loss of cellular identity during ageing".

Thank you for the suggestion – we have now toned down the title: ‘Inter-tissue convergence of gene expression during ageing suggests age-related loss of tissue and cellular identity’

Reviewer #3 (Recommendations for the authors):

Authors only analyze the gene sets that are correlated with age i.e. linearly go up or down with age. However, ageing trajectories for different genes can be very dynamic. This is the likely reason why only a few age-related changes were detected in some tissues. To assess the robustness of some of the key observations, authors should consider an orthogonal and unbiased clustering analysis similar to Figure 2 supplement to get the patterns for gene expression changes with age and to what extent DiCo is observed with different clusters.

We thank the reviewer for the comment. The reason we used Spearman’s correlation was to capture potential non-linear but monotonic changes. However, we agree with the reviewer that Spearman’s correlation cannot identify non-monotonic changes. As suggested, we now performed k-means clustering for each tissue to identify age-related expression patterns and then performed an enrichment test to assess if DiCo genes are enriched or depleted in any cluster. The results are now presented as supplementary figures (Figure 1—figure supplement 12-15). Notably, we did not observe a bias towards a cluster having a linear pattern with age or not in any of the tissues. For example, the cortex tissue displayed 15 clusters of expression changes, among which 5 are enriched and another 5 are depleted in DiCo (Figure 1—figure supplement 12). Two of the DiCo enriched clusters (2/5) do not have linear trajectories with age (clusters 12 and 14) whereas three of the DiCo depleted clusters (3/5) display linear patterns (clusters 6, 7 and 15). We have now presented this result in the relevant section (lines 276-277).

Some of the key observations have either low effect size or are not statistically significant. For example, inter-tissue expression similarity during ageing and development (lines 115-119), the expression reversal (lines 158-161), changes in coefficient of variation with age (lines 202-207) etc. Authors clearly mention these in the manuscript, but it raises concerns on the extent to which some of these observations might be random, or due to idiosyncrasies of RNA-seq itself. These caveats should also be discussed.

We now added a limitations paragraph in the Discussion, detailing both this aspect and other limitations raised by the other reviewers (lines 605-620).

In the global PCA based analysis for figure 1 and 2 combining both development and ageing data, the variance can be driven by one of the two states. For example, transcription divergence analysis in PCA in figure 1 based on development alone is not significant which suggests that the divergence of the developmental samples is driven in-fact by ageing. For PCA based analysis using separate PCAs for development and ageing will be more statistically sound.

We thank the reviewer for this suggestion. We now followed the reviewer’s suggestion and analysed age-related convergence using PCA performed only on ageing samples. In this way, we now analyse Euclidean distance in PC1-4 (i) using all samples together (across the lifetime – Figure 1a), (ii) using samples from developmental period only (Figure 1—figure supplement 3a-b), and (iii) using samples from ageing period only (Figure 1—figure supplement 3d-e). For the PCA analysis based on all the samples together, change in mean Euclidean distance with age among tissues was significant both during development and ageing (⍴_dev=0.99, p_dev=1.5x10^-5; during ageing ⍴_age=-0.87, p_age=0.0026). When we performed the PCA analysis separately for the two periods, we still observed a significant divergence among tissues during development (⍴=0.95, p=0.0008) and observed a marginally significant convergence during ageing (⍴=-0.64, p=0.0596). We hope that this analysis addresses the reviewer’s point. We also note that our earlier analysis had an error (we had used PC3-4 spaces to calculate the Euclidean distance for the mouse dataset, mistakenly excluding PC1-2), which we now corrected (Figure 1, Figure 1—figure supplement 3).

Are there any known gene categories that have stronger expression reversal trends? For example, do known transcription factors or pathways that regulate development show DiCo if analyzed separately? Or is this trend only seen at the genome-wide level and completely stochastic?

We thank the reviewer for the suggestion. We now used the MiRTarBase database for miRNAs and TRANSFAC for transcription factors and tested if any specific regulator is associated with DiCo pattern genes. We did not find any significant association (lines 353357). This may be consistent with DiCo being a transcriptome-wide and diffuse phenomenon caused by stochastic changes, as the reviewer indicates. It is also possible that not specific regulators but a change in their cooperative action influences the pattern, or we lack sufficient power with the current datasets to identify subtle effects of individual regulators. We now mention this also in the Discussion (line 601-603).

For functional enrichment analyses authors only mention some selected functions in the text but GO terms have a lot of redundancy and it is unclear if there is any emerging functional trends looking beyond the redundancy or selected examples. It will be useful to thoroughly analyze and summarize the enriched GO terms to assess for any functional patterns that connect functional decline with ageing and convergent expression or DiCo.

In the previous version, enriched GO terms were summarised according to their Jaccard similarities using the ‘emapplot’ function of the ‘enrichplot’ package. However, we agree with the reviewer that the enrichment network did not remove redundancy among the categories, and might have hidden unique and interesting functions.

We have now used another approach to summarise and interpret DiCo enriched categories in more detail. First, we clustered the categories based on the number of genes they share using hierarchical clustering and cut the tree into 25 clusters. For each cluster, we chose a representative category that has the highest mean-Jaccard similarity to the other categories in the same cluster (Dönertaş et al., 2021). Then, we calculated the mean expression changes (⍴ between expression and age) of the genes in the representative categories for each tissue, separately (Figure 4h). We believe that this new analysis provides better insight to our main findings such that tissue-related functions gained during development display opposing expression trends during ageing, i.e. reversal pattern, providing support to loss of tissuespecific expression and thus contributing to convergence during ageing. We have now updated the main text and methods (lines 343-351, 874-881).

Is there a significant overlap between mouse and human ageing related changes at the gene or pathway level?

We thank the reviewer for the suggestion. We now compare the correlations between ageing related expression changes and overlap between convergent genes during ageing (Figure 4 figure supplement 2). While mouse datasets, especially our data and Jonker et al., showed higher correlations among each other, the correlations between age-related expression changes in the human and mouse datasets were smaller but overall positive when we consider the same tissues in both datasets (average correlations between GTEx and mouse datasets are ⍴_Jonker=0.12, ⍴_{Our_dat}=0.13, ⍴_Schaum=0.01). Convergent genes in ageing also showed weak but significant overlaps across datasets. Human dataset showed significant overlap with Jonker et al., and Schaum et al., datasets but not with ours. Nevertheless, the overlap for these datasets was also small. Importantly, apart from the differences in species, sex, age distribution and the technology used to generate the data, the small overlap between datasets may stem from the fact that the external datasets lack developmental samples. Thus, we can only compare convergence during ageing but not the DiCo pattern, the main focus of our study. Since only 62% of the convergent genes during ageing are divergent in development in our dataset, we should emphasise that low overlap for convergence does not rule out overlap across DiCo genes. Similarly, we use divergent genes during development as the background for our functional enrichment tests to find functional associations of DiCo. Since we lack the developmental background we performed functional enrichment for only the convergent genes but did not find any significant association for the GTEx dataset.

Line 299 – linking the functional categories with cellular identity appears a bit hand-wavy. Stronger support would be needed to make this conclusion.

Thank you for the suggestion. We have now toned down that particular line (347-350), and also revised the manuscript, in general, to make sure we do not overstate our findings.

Some of the aspects need more methodological details:

For RNA-seq analysis, were all the samples normalized together, or were development and ageing samples normalized separately? It should be clarified in the methods.

We thank the reviewer for pointing out this ambiguity. We have combined all the samples together regardless of their age (development and ageing together) or tissue (n=63) and performed quantile normalisation on log2 transformed (after adding 1) FPKM values. As FPKM normalisation does not account for inter-sample variability, we adopted this approach to control for cross-library variability between different tissues and ages. We have rewritten the relevant section in the methods to explain the normalisation in more detail (line 666-685). The same approach is used for the VST-normalisation in this version and development and ageing periods are normalised together.

Were the mouse perfused? The details should be added in methods. If the mouse were not perfused, a significant shared proportion of gene expression can be due to blood components. Also, the possibility that some of the shared changes could be due to infiltration of immune and other cell types in tissues should be mentioned.

We thank the reviewer for this question. Indeed, infiltration of other cell types is a serious consideration in most bulk tissue analyses. Since we observe a similar trend between different cell types using scRNAseq data and that DiCo is particularly associated with tissue-specific genes and not with any bood or immune-related signature, we believe the effect should be minimal in our analysis. However, the new dataset we included in the analysis, Schaum et al., did use perfused mice and although the association between the loss of tissue identity and convergence was particularly strong, the transcriptome-wide trend did not support convergence. We now included the potential influence of infiltration of other cell types as a limitation of our study in the Discussion (line 606-612).

More details for single cell deconvolution should be added in methods.

We have now updated the deconvolution section in the methods to explain how we conducted the analysis in more detail (lines 1054-1064).

It is unclear to me how how the age related gene expression and PC correlation was performed. Was the loadings underlying PC axes used for correlation analysis? It should be elaborated in methods.

We performed PCA on individuals’ gene expression values and used these PC scores of the individuals (not the loadings) and their ages to perform the correlation analysis. We have now updated the relevant section in the methods (lines 693-696).

Some of the supplemental figures are missing.

We are sorry for this confusion, but having checked our submission we could not detect any missing files. In case this might be an issue with the submission portal, for the second round we will update the bioRxiv version of the article where the reviewers can check all files separately if needed.

Author details

1. Hamit Izgi

   Department of Biological Sciences, Middle East Technical University, Ankara, Turkey

   Contribution

   Conceptualization, Data curation, Formal analysis, Methodology, Software, Visualization, Writing – original draft, Writing – review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-4030-3132

2. Dingding Han

   CAS Key Laboratory of Computational Biology, CAS-MPG Partner Institute for Computational Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China

   Present address

   Department of Clinical Laboratory, Shanghai Children's Hospital, Shanghai Jiaotong University, Shanghai, China

   Contribution

   Investigation, Resources, Writing – review and editing

   Competing interests

   No competing interests declared

3. Ulas Isildak

   Department of Biological Sciences, Middle East Technical University, Ankara, Turkey

   Contribution

   Data curation, Formal analysis, Validation, Visualization, Writing – review and editing

   Competing interests

   No competing interests declared

4. Shuyun Huang

   CAS Key Laboratory of Computational Biology, CAS-MPG Partner Institute for Computational Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China

   Contribution

   Investigation, Writing – review and editing

   Competing interests

   No competing interests declared

5. Ece Kocabiyik

   Department of Biological Sciences, Middle East Technical University, Ankara, Turkey

   Contribution

   Data curation, Formal analysis, Writing – review and editing

   Competing interests

   No competing interests declared

6. Philipp Khaitovich

   Center for Neurobiology and Brain Restoration, Skolkovo Institute of Science and Technology, Moscow, Russian Federation

   Contribution

   Conceptualization, Project administration, Resources, Supervision, Writing – review and editing

   For correspondence

   p.khaitovich@skoltech.ru

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-4305-0054

7. Mehmet Somel

   Department of Biological Sciences, Middle East Technical University, Ankara, Turkey

   Contribution

   Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review and editing

   For correspondence

   somel.mehmet@googlemail.com

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-3138-1307

8. Handan Melike Dönertaş

   1. European Molecular Biology Laboratory, European Bioinformatics Institute EMBL-EBI, Wellcome Trust Genome Campus, Cambridge, United Kingdom
   2. Leibniz Institute on Aging - Fritz Lipmann Institute (FLI), Jena, Germany

   Contribution

   Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing

   For correspondence

   melike.donertas@leibniz-fli.de

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-9788-6535

• 3,650

  Page views

• 465

  Downloads

• 11

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.