Inter-tissue convergence of gene expression and loss of cellular identity during ageing

124/150 words) Gene expression tends to reverse its developmental trajectory direction in ageing, a phenomenon previously linked to cellular identity loss. Our analysis of cortex, lung, liver and muscle transcriptomes of 16 mice, covering development and ageing intervals, revealed widespread but tissue-specific ageing-associated expression reversals. Cumulatively, these reversals create a unique phenomenon: mammalian tissue transcriptomes diverge from each other during postnatal development, but during ageing, they tend to converge towards similar expression levels, a process we term divergence convergence, or DiCo. We confirmed the DiCo pattern using data from independent mouse and human datasets. Further, using publicly available single-cell transcriptome data, we showed that DiCo could be driven both by alterations in tissue cell type composition and also by cell-autonomous expression changes within particular cell types.


Introduction
Development and ageing in multicellular organisms are highly intertwined processes. On the one hand, certain ageing-related phenotypes, such as presbyopia, have been suggested to be caused by the continuation of developmental processes into adulthood (Blagosklonny 2013;de Magalhães and Church 2005). Such "runaway development" phenomena may be due to insufficient natural selection pressure to optimise regulation after sexual reproduction starts. Conversely, 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 relative reduction in neuronal proportions (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 toward 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 needs to be explored.
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 hyperfunction of developmental genes. However, such work has been limited by the scarcity of studies that comprise 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.

Results
We generated bulk RNA-seq data from 63 samples covering the cerebral cortex (which we refer as cortex), liver, lung, and skeletal muscle (which we refer as muscle) of 16 male C57BL/6J mice, aged between 2 to 904 days of postnatal age (see Methods). As mice reach sexual maturity by around two 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, M. Currer, and Harrison 2007)) (n=9) as the ageing series (Figure 1figure supplement 1). The final dataset contained n=15,063 protein-coding genes expressed in at least 25% of the 63 samples (one 904 days 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 -16 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 correlation coefficient ⍴=[0.88, 0.99], nominal p<0.01 for each test; Figure 1b). Also, after the tissue effect was removed by standardisation, the PCA showed a strong influence of age on the first two PCs, which explain 31% of the variance in total (Figure 1-figure supplement 2).
In the PC3-PC4 space of Figure 1a, we further observed a trend of higher similarity among juvenile samples than among adult samples across tissues (change in mean Euclidean distance among tissues with age during development ⍴dev=0.90, pdev=0.006; during ageing ⍴age=-0.30, page=0.429, 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 principal components analysis (PCA) was performed with developmental samples only but was not significant (days 2 to 61; Figure 1-figure supplement 3).
Tissues involve common gene expression changes with age. We next characterised the agerelated changes in gene expression shared across tissues by studying i) the trends at the transcriptome level and ii) the significant changes at the gene level. First, we studied similarities in overall trends in age-related gene expression changes using the Spearman correlation coefficient (⍴) between expression levels and age, for each gene, in each tissue, separately for the developmental and ageing periods (see Methods, tissue-specific age-related gene expression changes and functional enrichment test results are available as Supplementary File 1). We then examined similarities across tissues during development and ageing by comparing the gene-wise expression-age correlation coefficients (Figure 1c). Considering the whole transcriptome without a significance cutoff, we found a weak correlation of age-related expression changes in tissue pairs, both during development (⍴=[0.17, 0.39], permutation test p<0.05 for all the pairs, Figure 1-source data 1), and ageing (⍴= [0.23, 0.33], permutation test p<0.05 in 4/6 pairs, Figure 1-source data 1). The difference between inter-tissue similarities within the development and those within ageing was not significant (Wilcoxon signed-rank test, p=0.31). Moreover, the number of genes with the same direction of change across four tissues was consistently more than expected by chance (p<0.05), except for genes upregulated in ageing (Figures 1e, Figure 1-figure supplement 4). This attests to overall similarities across tissues both within the development and within ageing datasets, albeit of modest magnitudes.
In the second approach, we focused on genes with a significant age-related expression change, identified separately during development and ageing (using Spearman correlation and false-discovery rate (FDR)<0.1, Figure 1d). We found that the developmental period was accompanied by a large number of significant changes (n= [1,941,6,151], 13-41% across tissues), with the most manifest changes detected in the cortex. The genes displaying significant developmental change across all four tissues also showed significant overlap (Figure 1-figure supplement 5a, Figure 1-figure supplement 6; permutation test: pshared_up=0.027, pshared_down<0.001, false-positive rate < 20%). Using the Gene Ontology (GO), we found no functional enrichment for these shared developmentally upregulated genes (Benjamini Yekutieli (BY)-corrected p<0.1). In contrast, shared developmentally down-regulated genes were enriched in functions such as cell cycle and cell division (BY-corrected p<0.1; Supplementary File 2). Contrary to widespread expression change during development (13-41%), the proportion of genes undergoing significant expression change during ageing was between 0.0001-15% (Figure 1d). This contrast between postnatal development and ageing was also observed in previous work on the primate brain (Somel et al. 2010;Işıldak et al. 2020). In terms of the number of genes with the significant ageing-related change, the most substantial effect we found was in the lung (n=2,319), while close to no genes showed a statistically significant change in the muscle (n=2), a tissue previously noted for displaying a weak ageing transcriptome signature across multiple datasets (Turan et al. 2019). Not unexpectedly, we found no common significant ageing-related genes across tissues (Figure 1-figure supplement 5a). Considering the similarity between the ageing and development datasets (Figure 1c) and the similar samples sizes in development (n=7) and ageing periods (n=9), the lack of overlap in significant genes in ageing might be due to low signal-to-noise ratios in the ageing transcriptome, as ageing-related changes are subtler compared to those in development ( Figure 1-figure supplement 5b).
Gene expression reversal is a common phenomenon in multiple tissues. We then turned to investigate the prevalence of the reversal phenomenon (i.e. an opposite direction of change during development and ageing) across the four tissues. We first compared age-related expression change between development and ageing periods in the same tissue (Figure 1c). This revealed weak negative correlation trends in liver and muscle (though not in the lung and cortex), i.e. genes up-or down-regulated during development tended to be down-or up-regulated during ageing, respectively. These reversal trends were comparable when the analysis was repeated with the genes showing relatively high levels of age-related expression change (|⍴|>0.6; Figure 1-figure supplement 7). We further studied the reversal phenomenon by classifying each gene expressed per tissue into those showing up-or down-regulation (irrespective of statistical significance) during development and during ageing, thus summarising trends of continuous change versus reversal in each tissue. In line with the above results, as well as earlier observations in the brain, kidney, and liver (Dönertaş et al. 2017;Anisimova et al. 2020), we found that ~50% (43-59%) of expressed genes showed reversal trends (Figure 1f), although these proportions were not significantly more than randomly expected in permutation tests (Figure 1-figure supplement 8, see Methods). Overall, we conclude that the expression trajectories of the genes do not necessarily continue linearly into the ageing period.
Pathways related to development, metabolism and inflammation are associated with the reversal pattern. We then asked whether genes displaying reversal patterns in each tissue may be enriched in common functional categories. Our earlier study focusing on different brain regions had revealed that up-down genes, i.e. genes showing developmental up-regulation followed by downregulation during ageing, were enriched in tissue-specific pathways, such as neuronal functions (Dönertaş et al. 2017). Analysing up-down genes compared to all genes up-regulated during development, we also found significant enrichment (BY corrected p-value<0.1) in functions such as "synaptic signaling" in the cortex, as well as "tube development" and "tissue morphogenesis" in the lung, "protein catabolic process" in the liver and "cellular respiration" pathways in the muscle (Supplementary File 3). Meanwhile, down-up genes (down-regulation during development followed by up-regulation during ageing) showed significant enrichment in functions such as "wound healing", and "peptide metabolic process" in the cortex, "translation" and "nucleotide metabolic process" in the lung, and "inflammatory response" in the liver (Supplementary File 3).
Genes showing reversal pattern are not shared among tissues. As tissues displayed modest correlations in development-and ageing-related expression change trends (Figures 1c, Figure 1figure supplement 7), and as we had previously observed that distinct brain regions show similarities in their reversal patterns (i.e. the same genes showing the same reversal type), different tissues might also be expected to show similarities in their reversal patterns. Interestingly, we found no overlap 6 between gene sets with the reversal pattern (up-down or down-up genes) across tissues, relative to random expectation (pup-down>0.07, pdown-up=0.53, Figure 1-figure supplement 9). Such a lack of overlap might be explained if genes showing reversal patterns in each tissue tend to be tissuespecific. It would also be consistent with the notion that reversals involve loss of cellular identities gained in development, during which tissue transcriptomes appear to diverge from each other (Figures 1a, Figure 1-figure supplement 3) (Cardoso-Moreira et al. 2019). This result led us to ask whether, in accordance with the reversal phenomenon, inter-tissue transcriptome divergence may be followed by increasing inter-tissue similarity, or convergence, during ageing.
Inter-tissue divergence during development and convergence during ageing. We studied the inter-tissue divergence/convergence question using two approaches. In the first, we studied how expression variation among tissues changes with age. 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 correlation between CoV and age, separately for development and ageing periods (List of genes are given in suggesting the lack of significant result with all tissues might be due to decreased power. 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,2). Figure 2b exemplifies this pattern of increasing and then decreasing CoV through lifetime for the gene displaying the strongest such 7 signal. We next studied the transition points between divergence and convergence by clustering genes showing divergence-convergence pattern (without using a significance cutoff, n=4,802) based on their CoV change with age ( Figure 2-figure supplement 3). Notably, the cluster with slightly delayed divergence (Cluster 1) was associated with metabolic and respiration-related processes (BYcorrected p<0.1), and the cluster with a relatively delayed convergence (Cluster 5) was enriched in categories related to vascular development (BY-corrected p<0.1) (Supplementary File 4). To assess the contribution of tissues to divergence-convergence pattern, we clustered the genes (n=4,802) based on their expression change with age ( Figure 2-figure supplement 4). Not surprisingly, the clusters with relatively higher expression levels of a tissue (e.g. muscle in Cluster 9) were enriched in functional categories (BY-corrected p<0.1) related to that tissue (e.g. muscle cell development)

(Supplementary File 5).
We then studied this pattern at the single-gene level. We tested each gene for significant CoV change in their expression levels (i.e. divergence or convergence) in development and ageing (Spearman correlation test with FDR<0.1). We found that the ratio of divergent and convergent genes differed significantly between development (70% divergence among 2,581 significant genes) and ageing (68% To our knowledge, inter-tissue convergence during ageing is a novel phenomenon. In order to verify this trend, we used a second approach to test it in our mouse dataset. Inter-tissue convergence during ageing is a reproducible trend in another mouse dataset and in humans. We investigated inter-tissue convergence during ageing in another mouse ageing transcriptome dataset that comprises 5 tissue samples (brain, lung, liver, kidney, spleen) collected from the same mice (n=18) (Jonker et al. 2013). We confirmed the decrease in mean CoV during ageing (⍴ = -0.58, p = 0.01) as well as the increased pairwise tissue correlations in 7/10 tissue pairs (Figure 2-figure supplement 8).
We next investigated inter-tissue convergence in the GTEx dataset, which includes transcriptome data from humans covering a large set of tissues during the ageing period (GTEx Consortium et al. 2017). We first confirmed the results using 47 individuals with samples in the same four tissues as those studied in mice. We first noticed a tendency for an age-related increased similarity between tissue samples of the same individual using the mean Euclidean distance between PCs ( ⍴ = -0.23, p = 0.12, Figure 2-figure supplement 9). We then calculated mean and median CoV values (⍴meanCoV = -0.12, pmeanCoV = 0.41; ⍴medCoV = -0.18, pmedCoV = 0.22) and pairwise inter-tissue correlations (5/6 comparisons positive, none was significant); which also followed a trend compatible with intertissue convergence during ageing, albeit not significant (Figure 2-figure supplement 10). We further tested inter-tissue convergence using 10 tissues in GTEx, covering 35 samples from adipose, tibial artery, cerebellum, lung, skeletal muscle, tibial nerve, pituitary, sun-exposed skin, thyroid, and whole blood (see Methods). Similar to our previous results, a tendency for inter-tissue convergence during ageing was observed using through PCA (⍴ = -0.26, p = 0.13, These results suggest that inter-tissue convergence in ageing is a weak but widespread phenomenon, observed in multiple tissues in mice and possibly also in humans. Overall, while mouse and human tissues display divergence in development (Figures 1a, 2a, (Cardoso-Moreira et al. 2019)) this is followed by convergence in ageing (Figures 2a, Figure 2-figure supplement 1-12). This we refer to as the divergence-convergence (DiCo) pattern of postnatal lifetime.

The divergence-convergence (DiCo) pattern indicates loss of tissue-specificity during ageing.
Potential explanations of DiCo pattern involve two scenarios consistent with the age-related loss of cellular identity: i) decreased expression of tissue-specific genes in their native tissues, or ii) nonspecific 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: 1,175, lung: 839, liver: 986, muscle: 766 genes). We noted that tissue-specific genes We then tested our initial prediction that the DiCo pattern is related to tissue-specific genes losing their expression in their native tissue during ageing, or gaining expression in non-native tissues. We first tested the loss of tissue-specific expression during ageing by considering all tissue-specific genes and found a positive odds ratio (OR = 5.50, Fisher's exact test p=2.1x10 -129 , Figure 4a). The same analysis with only the DiCo genes yielded a much stronger association (OR=74.81, Fisher's exact test p=5.9x10 -203 , Figure 4b). This suggests that loss of tissue-specific expression is observed across the transcriptome, with a particularly strong association among DiCo genes.
We then asked whether genes displaying the DiCo pattern may be related to specific functional pathways. Using GO, we searched for functional enrichment among convergent genes during ageing, using developmentally divergent genes as the background (Methods). We found enrichment for 237 GO Biological Process (BP) categories for the DiCo pattern (Kolmogorov-Smirnov (KS) Test, BYcorrected p<0.1, Figure 4-source data 1). In addition to mitochondria, energy metabolism and other metabolic pathways, we also found categories related to tissue development and differentiation Cellular composition changes and cell-autonomous expression changes can both explain the divergence-convergence pattern. Ageing-related transcriptome changes observed using bulk tissue samples may be explained by temporal changes in cell type proportions within tissues, by cellautonomous 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 singlecell 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 (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=[4,007,4,106]) and non-DiCo genes (n=[8,485, 8,743]). 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 Hence, these results indicate that the observed cellular composition changes partly explain but are not exclusive to genes displaying the DiCo pattern. Next, we investigated the possible role of cell-autonomous changes in the DiCo pattern. Cellautonomous changes could contribute to inter-tissue convergence during ageing in two ways. First, similar or corresponding cell types present in different tissues, such as immune cells, might become more similar with age in their expression profiles. 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, become more similar, i.e. converge, with age. To test these scenarios, we first ordered the pairwise correlations between cell types in different tissues to determine the most similar and dissimilar cell types across tissues (Methods). Then, we studied how these similarities (i.e. pairwise correlations) change with age ( Figure 5b). Intriguingly, we found that similar cell types (i.e. those with the highest correlations) among tissues 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 re-sampling all genes, p>0.1; and 47/54 [87%] increase in correlation among the most distinct cell types, permutation test, p>0.1). These results suggest that cell-autonomous changes also contribute to inter-tissue convergence during ageing and provide evidence for age-related cellular identity loss.
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 tissue from 3m-to 24m-old mice (Figure 5-figure supplement 6).

Discussion
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 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), 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). Parallel expression changes among tissues during ageing may be related to common responses to damage and inflammation and processes such as reduced energy metabolism, as observed earlier (Zahn et al. 2007).
Second, we verify the generality of the reversal pattern, i.e. 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. We thus propose that reversal is a general phenomenon in mammalian ageing.
Two observations are notable. One is that reversal-displaying genes, especially those displaying the up-down pattern in each tissue, can be associated with tissue-specialization-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. These imply that the parallel changes among tissues observed during development and ageing, described earlier, do not include reversals. Therefore, we hypothesised that reversals might instead be linked to inter-tissue divergence and tissue specialisation during development, and possibly convergence during ageing.
We indeed observed that the tissue-specific genes show an up-down reversal pattern in their native tissue. 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, albeit weak, towards inter-tissue convergence during ageing. Studying two additional datasets from mice and humans, we observed that convergence during ageing is common, although possibly not ubiquitous.
Importantly, these results reliably show that developmental inter-tissue divergence does not continue into ageing. We predict that future work on ageing with datasets, including multiple distinct tissues and larger sample sizes may establish the prevalence of the DiCo phenomenon.
Finally, we report three interesting observations on DiCo. We determine that tissue-specific genes tend to be down-regulated in the tissues that they belong to during ageing, while non-tissue-specific 13 genes are up-regulated. 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 (Tabula Muris Consortium 2020). Accordingly, the DiCo phenomenon can 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 growing expression similarity between normally distinct cell types in different tissues during ageing, but decreased similarity between normally similar cell types.
Cell-autonomous expression changes, therefore, may also contribute to the divergence-convergence phenomenon.
These results support a model where expression reversals and the newly discovered divergenceconvergence phenomenon are associated with loss of specialisation, at the tissue level and possibly also at the cellular level.

Sample Collection
We collected time-series 16 male C57BL/6J mice samples representing the whole lifespan of Mus musculus, comprising both postnatal development and the ageing period from four different tissues; cerebral cortex, liver, lung and skeletal muscle. One 904 days-old mouse had no cortex tissue sample, and thus excluded from the analysis. As a result, we have generated 63 RNA-seq libraries in total.

Separation of development and ageing periods:
In order to compare the gene expression changes during postnatal development and ageing we studied the samples before sexual maturation (covering 2 to 61 days of age , n=7) as postnatal development period, and samples covering 93 to 904 days (n=9 in all except in cortex n=8) as the ageing period.

RNA-Seq Library Preparation
RNA sequencing was performed as previously described (Liu et al. 2016) with slight modifications.
Briefly, total RNA was extracted using Trizol reagent (Invitrogen). For sequencing library construction, all samples were randomized first and prepared with TruSeq RNA Sample Preparation Kit (Illumina) according to the manufacturer's instruction. Libraries were then sequenced on the Illumina Hi-seq 4000 system in three lanes within one flow-cell, using the 150-bp paired-end module.

RNA-Seq Data Preprocessing
The quality assessment of the raw RNA-seq data was done using FastQC (Andrews 2010) and the adapters were removed using Trimmomatic (Bolger, Lohse, and Usadel 2014). The low-quality reads were filtered using the parameters: PE ILLUMINACLIP: TruSeq3-PE-2.fa:2:30:1:0:8:true, SLIDINGWINDOW:4:15, MINLEN:25. Remaining high-quality reads were aligned to the reference genome GRCm38 using STAR-2.5.3 (Dobin et al. 2013) with parameters: --sjdbOverhang 99 --outSAMattrIHstart 0 --outSAMstrandfield intronMotif --sjdbGTFfile GRCm38.gtf. Percentage of uniquely mapped reads in libraries ranged from 80 to 93%. We used cufflinks (Trapnell et al. 2010) to generate read counts for uniquely aligned reads (samtools -q 255 filter) and calculated expression levels as fragment per kilobase million (FPKM). There were 50 duplicated genes, and the sum of their FPKM values were used. In total we had 51,608 genes. We restricted the whole analysis to only protein-coding genes obtained by the 'biotype' feature of the biomaRt library (Durinck et al. 2009). We also excluded genes which were not detected (zero FPKM) in 25% or more of the samples (at least 15 of 63), resulting in 15,063 protein-coding genes in total. Using all the samples together, FPKM values were log2 transformed (after adding 1) and quantile normalized with 'normalize.quantiles' function from 'preprocessCore' library (Bolstad 2020) to be used in downstream analysis.

Principal component analysis:
We estimated the source of variation in the data using Principal component analysis (PCA) on the scaled expression matrix with 'prcomp' function in the R base. The first four components, from PC1 to PC4 in order, explained 31%, 20%, 17% and 8% of the total variance with a clear separation of tissues in PC1 and PC2 and a strong age effect in PC4. We performed ANOVA between each PC1-3 and tissue labels to confirm tissue differences. The magnitude of age effect on PCA analysis was measured with Spearman correlation between individual age and each principal component. We further calculated euclidean distance in pairwise manner among tissues of each individual in PC3-4 space constructed using only the developmental samples. Then, we tested the effect of age on mean euclidean distance among tissues using Spearman correlation.

Age-related gene expression change
To identify genes showing age-related expression change in each tissue, we used Spearman correlation between individual age and expression level, separately for development and ageing periods. To capture potential non-linear but monotonic changes in expression, we chose the nonparametric Spearman correlation test for both periods. Significance of age-related genes was assessed with FDR<0.1 cutoff, calculated with the Benjamini-Hochberg (BH) procedure (Benjamini and Hochberg 1995) using the 'p.adjust' function in the R base library.

Functional associations:
We tested the functional associations of age-related gene expression change in separate tissues for each period separately, following gene set over-representation analysis (GORA) procedure with Gene Ontology (GO) (Ashburner et al. 2000) Biological Process (BP) categories using the 'topGO' package (Alexa and Rahnenfuhrer 2019). We applied the classical algorithm and performed Fisher's Exact Test on categories that satisfy the criteria of a minimum 10 and maximum 500 number of genes. We used the whole set of expressed genes (n=15,063) as the background. P-values were corrected for multiple testing using Benjamini-Yekutieli (BY) (Benjamini and Yekutieli 2001) procedure, which does not assume independent tests and thus is suitable for multiple testing correction for GO enrichment with potentially overlapping gene sets. Categories with BY corrected p-value<0.1 were considered as significant.

Correlation between age-related gene expression changes in different tissues
We calculated Spearman correlation coefficients between age-related gene expression change values in each tissue pair (⍴ from the correlation between gene expression levels and age) (Figure 1c). In order to test the statistical significance of the correlations we used a permutation scheme as the expression levels across tissues are not independent but of the same mice. In order to account for the dependence, the individual ages were permuted but the permuted values kept the same across tissues (similar to permutation tests applied in (Dönertaş et al. 2017;Işıldak et al. 2020;Dönertaş et al. 2018)). Specifically, we randomized the individual ages for 1000 times, in development and ageing separately, using the 'sample' function in R, while keeping the permuted age labels the same for individuals across tissues. We calculated the age-related gene expression changes with permuted ages, thus simulating the null distribution with no age effect.
We then calculated the Spearman correlation between the age-related gene expression levels from the permutations across tissues and we assigned the p value by calculating the proportion of permuted calculations with a more extreme correlation. FPR or false positive rate was calculated as the median value of expectations divided by the observed value (Figure 1-source data 1).

Shared gene expression changes across tissues:
We summarized the number of shared age-related genes among tissues in both development and ageing and for up-and down-regulated genes separately, using FDR<0.1 (Figure 1-figure   supplement 5). For each gene, we counted the number of tissues with the same direction of expression change with age. We also calculated the same overlap among tissues with all genes without using any significance cutoff (Figure 1e).

Permutation test:
We again used a permutation scheme to assess the significance of shared age-related genes to account for the dependence between tissues. We tested the significance of shared up-and down-regulated genes, selected with or without FDR cutoff, and development and ageing periods separately. We used the permuted age-related expression changes that are described before for the test of correlations. To test the significance of the overlap of significantly up-or down-regulated genes (FDR<0.1), we ranked the permuted age-related expression changes (⍴) for each tissue in each period separately, and chose the highest N (to test the up-regulation), or lowest N (to test the down-regulation) number of genes, where N is the number of significantly up-or downregulated genes in a given tissue (FDR<0.1). We then calculated the number of overlaps across tissues and thus generated a distribution of shared expression changes when there is no age effect. We calculated the p value as the number of the permuted values with a higher number of overlap than the observed value divided by the number of permutations (n=1,000). FPR was calculated as the median number of shared expression changes in permutations divided by the observed value.

Functional Associations:
We tested the functional associations of shared expression change trends among tissues in each period, separately, following the GORA procedure using the same criteria and algorithms explained in the previous section. To test shared up-regulated (n=45) or down-regulated genes (n=138) in development, we chose all significant age-related genes across tissues (n=10,305) in the development period. Since we could not identify any shared ageing-related genes across tissues (Figure 1-figure supplement 5), we did not perform a functional test for the ageing period.

Analysis of gene expression reversals
We compared the direction of gene expression change during development and ageing to identify

Permutation test:
To test the significance of reversal proportions, we kept the developmental changes constant and randomly permuted the individual ages only in the ageing period (as described earlier). Among developmental up-regulated genes, we calculated an UD% from each permutation, simulating a null distribution for UD reversal. We applied the same principle for the DU genes. Thus, we created a null distribution with the expected reversal ratios and tested the significance of observed values for each tissue separately.

Functional associations:
We used the GORA procedure as described earlier but kept the developmental changes constant in the background. More specifically, we tested the functional enrichment of UD reversal genes against UU genes, and DU genes against DD genes to keep the developmental changes fixed and test the reversal pattern not developmental associations.

Overlap of reversal genes -permutation test
We tested the significance of overlap using the same permutation scheme described above.
Specifically, among developmental up-(or down-) regulated genes shared among tissues (nup=2,255, ndown=2,209), we constructed the null distribution by calculating the ratio of UD vs UD+UU (or DU vs DU+DD) genes shared among tissues in ageing permutations (Figure 1-figure supplement 8).

Tissue convergence and divergence calculations using coefficient of variation (CoV)
For each individual mouse and gene (n=15,063), we calculated a coefficient of variation (CoV) with normalized expression levels from the four tissues by dividing the standard deviation by the mean. We studied the (inter-tissue) expression-variation change with age in development and ageing periods separately, using two approaches.

Mean/median CoV across genes:
We assessed the global variation among the tissues of each individual mouse by calculating the mean (or median) CoV of genes and then using Spearman correlation test between mean-CoV (or median-CoV) and individual age.

CoV at the gene level
In the second approach, we tested each commonly expressed gene using Spearman correlation test between the CoV value of a gene and age. P-values were corrected for multiple testing, using the 'BH' procedure and the significance was assigned to corrected p-value<0.1. The genes showing positive correlation between CoV and age were called divergent, and the ones showing negative correlation between CoV and age were called convergent (Figure 2d). Genes that show divergent pattern during development and convergent pattern in ageing (without using a significance level) were called divergent-convergent (DiCo) genes (n=4,802).

Permutation Test
To test the significance of DiCo genes (n=4,802), we kept the developmental divergent genes constant (n=9,058, without a significance cutoff) and randomly permuted the individual ages only in the ageing period (as described earlier). Among developmental divergent genes, we calculated DiCo % for each permutation, simulating a null distribution for DiCo pattern (Figure 2-figure supplement   13).

DiCo Genes Clustering
We used the k-means algorithm to cluster DiCo genes according to their CoV or expression changes with age, separately. To find the optimum number of clusters for both procedures, we applied gap statistics using the 'clusGap' function in the 'cluster' package with 500 simulations (Tibshirani, Walther, and Hastie 2001). We used the 'kmeans' function in base R with 'iter.max=20' and 'nstart=50' parameters to cluster CoV values or expression levels which were standardized to mean=1 and sd=0 across genes.

Functional association:
To test the functional associations of the genes showing divergent-convergent pattern among tissues, we performed GSEA using GO BPs. We retrieved developmental divergent genes (CoV-age ⍴>0, n=9,058) and multiplied these correlation values with the ones calculated in the ageing period. Therefore, the genes with a negative value represent a divergent-convergent pattern quantitatively, while the ones with a positive product represent a divergent pattern throughout lifetime. We then ranked the genes according to the calculated product values and sought enrichment for the upper and lower tail of the distribution using the Kolmogorov-Smirnov (KS) test implemented in the 'clusterProfiler' package (Yu et al. 2012). 'gseGO' function was used with parameters: nPerm=1000, minGSSize=10, maxGSSize=500 and pValueCutoff=1. Therefore, the enriched categories for the genes in the lower tail of the distribution would represent DiCo enrichment.

20
We further sought enrichment among DiCo genes that were clustered with the k-means algorithm in the previous step for both CoV and expression clusters, separately. Genes in each cluster were tested among all DiCo genes using the same GORA procedure as described before.

Jackknife to test the Di/Co ratio between dev and ageing
We tested the significance of divergent/convergent gene ratios using a jackknife resampling procedure in development and ageing periods, separately. Leaving out an individual in each iteration, we re-calculated the number of significant divergent and convergent genes and their ratios, denoted as si. Then, we constructed a 95% confidence interval using the following formula: var is the variance of divergent/convergent ratios (si) in leave-one-out resamples and s is the mean of si values.

Pairwise tissue divergence-convergence test
In order to further verify the inter-tissue divergent-convergent pattern that we observed between development and ageing periods, we used another approach based on expression correlations among tissues. We calculated pairwise Spearman correlation coefficients among tissues of the same individual mouse using all commonly expressed genes. For each tissue pair, we tested the effect of age on expression-correlations using Spearman correlation test in development and ageing periods, separately. We then calculated the mean (or median) of all six pairwise tissue correlations for each individual mouse (Figure 2-figure supplement 6). The same procedure was applied to other imported datasets (Figure 2-figure supplement 8, Figure 2-figure supplement 10).

Determination of tissue-specific genes
To identify which tissue(s) contribute to the reversal pattern, we assigned each gene to a tissue to identify tissue-specific expression patterns. First, we calculated an effect size (ES) between the expression of a gene in a tissue versus other three tissues using the development samples only, and repeated this procedure for all tissues. Hence, we obtained ES for each commonly expressed gene in each tissue. ES was calculated using the 'Cohen's d' formula defined as the difference between the two means divided by the pooled standard deviation. We then assigned each gene to a tissue in which the gene has the highest ES. Finally, we retrieved only the fourth quartile (>Q3) of genes assigned to a tissue to define tissue-specific expression. Using this approach, we found 3,766 tissuespecific genes in total (cortex: 1,175, lung: 839, liver: 986, muscle: 766 genes).

Enrichment test with the direction of age-related change
We tested the association between tissue-specificity and the age-related change in expression during ageing using Fisher's exact test. Specifically, we constructed the contingency table with two groups; the first group defines the direction of maximum expression change of each tissue specific gene (being either positive or negative), the second group defines whether this maximum expression change a tissue specific gene occurs in its native tissue or not (either yes or no). Hence, a test result where a positive OR suggest either the expression of genes decrease the most in their native tissue or increase the most in non-native tissues.

Enrichment in DiCo pattern
We tested the association between tissue-specificity and DiCo pattern using Fisher's exact test, calculating the enrichment of tissue-specific genes within DiCo genes.

Jonker
We downloaded the raw data from the GEO database with GSE34378 accession number (Jonker et al. 2013) and followed the same analysis pipeline described above using all the samples from 5 tissues ("Brain -Cortex", "Lung", "Liver", "Kidney", "Spleen") of 18 mouse comprising 90 samples in total. This dataset represents ageing period of the mouse ranging from 90 to 900 days. Using the oligo package (Carvalho and Irizarry 2010), we retrieved the expression matrices and removed the probesets that are annotated to more than one gene. We confined the analysis to only the proteincoding genes expressed in at least 25% of all samples. Resulting 16,647 genes were log2 transformed (after adding 1) and quantile normalized using the preprocessCore library (Bolstad 2020) across all samples. Downstream analysis was the same as described above.

GTEx
We downloaded GTEx v8 data (GTEx Consortium et al. 2017) from the data portal and repeated the analysis in human tissues. We first confirmed our results in the same 4 tissues ("Brain -Cortex", "Lung", "Liver", "Muscle -Skeletal") and then expanded the analysis to 10 tissues ("Adipose -Subcutaneous", "Artery -Tibial", "Brain -Cerebellum", "Lung", "Muscle -Skeletal", "Nerve -Tibial", "Pituitary", "Skin -Sun Exposed (Lower leg)", "Thyroid", "Whole Blood"). In order to choose which tissues to analyse, we first choose the minor tissues with the highest number of samples for each major tissue, which prevents the representation of the same tissue multiple times. We then performed hierarchical clustering of tissues based on the presence of samples from the same individuals ( Figure   2-figure supplement 14) and cut the tree into 3 clusters based on visual inspection. We selected the cluster with the highest number of overlapping individuals to analyse. The same procedure was followed for both 4-and 10-tissue analyses. In particular, we restricted the analysis to the individuals with samples in all tissues analysed and with a death circumstance of 1 (violent and fast deaths due to an accident) and 2 (fast death of natural causes) on the Hardy Scale (n = 47 for 4 tissue, n = 35 for 10 tissue). We removed duplicated genes from the analysis. Similar to our analysis with the mice data, we used only the protein-coding genes that are expressed in at least 25% of all samples, totalling 16,211 for 4 tissues and 16,305 for 10 tissues. The TPM values obtained from the GTEx data portal are log2 transformed (after adding 1), and quantile normalised using the preprocessCore library (Bolstad 2020) in R. Downstream analysis was the same as other datasets.

Preprocessing
We used Tabula Muris Senis dataset (Tabula Muris Consortium 2020) for scRNA-seq analysis as it is the only dataset to our knowledge that includes time-series samples covering old age, and the tissues present in our dataset. Seurat processed FACS data of the tissues: lung, liver, skeletal muscle and non-myeloid brain, were downloaded from figshare database (Pisco 2020). Seurat package (Stuart et al. 2019) was used to retrieve the expression matrix of the cells that are annotated to cell types in the original article. Each tissue contains samples from three time points; 90 (3m), 540 (18m) and 720 (24m) days old mice, totaling 14 samples each in lung, liver and brain and 9 samples in liver. We excluded cell types with less then 15 cells among all samples, and excluded genes if the expression level is 0 for all cells at a given age. This resulted in a median number of 99-382 cells assigned to cell types, 6-24 cell types and 16,951-22,122 genes across tissues. Using 3 month old mice, we calculated cell type specific expressions in each tissue such that we first calculated the mean expression levels among cells of an individual mouse for each cell type and then calculated the mean among individuals to obtain an average expression value for each cell type. Uniprot gene symbols were converted to Ensembl gene ids using biomaRt R package (Durinck et al. 2009).

Deconvolution
We used cell type specific expression profiles of 3-month-old mice to estimate relative contributions of cell types to the transcriptome profiles of tissues in our mouse dataset. We used a linear-regression based deconvolution method for each tissue using three genesets; all genes (n= [12,492,12,849]), DiCo genes (n=[4,007,4,106]) and non-DiCo genes (n=[8,485, 8,743]). Regression coefficients were used as relative contribution of cell types. We then tested the effect of age on cell type contributions using Spearman correlation test in development and ageing.

Cell type similarity & their change with age
To investigate the contribution of cell autonomous changes to inter-tissue convergence in ageing, we calculated pairwise cell type expression correlations among tissues and studied how these correlations change with age. Based on pairwise correlations in the 3-months age group, we identified the maximally and minimally correlated cell type pairs among tissues. Specifically, for each cell type in a given tissue (e.g. 10 cell types in the liver), we chose the minimally correlated cell types from the other three tissues (e.g. choose the minimally correlated cell type among 15 cell types in the cortex, choose the minimally correlated one among 24 cell types in the lung, and the same procedure for 6 cell types in the muscle) and repeated this procedure for the other three tissues, resulting in 54 cell type pairs. Then, we calculated Spearman correlation coefficients between age and minimally correlated cell type pairs identified in the 3-months age group. Likewise, we repeated the same 24 analysis for the maximally correlated cell type pairs among tissues.

Permutation Tests:
To test whether DiCo genes are significantly more associated with cell type proportion changes than non-DiCo genes, we performed a permutation test based on re-sampling procedure. For each tissue, we took random samples among all genes (n=[12,492, 12,849]) with size N (N being the number of DiCo genes in that tissue), and repeated deconvolution as explained above. By calculating cell type proportion changes with age for each random sample repeated 1000 times, we created the null distribution for each cell type. Then, we calculated the p values as the number of random samples having the same or higher cell type proportion change values divided by the observed value (cell type proportion changes with DiCo genes).
We applied a similar permutation scheme explained above to test cell type similarity change differences between DiCo and non-DiCo genes. For each random sample of non-DiCo genes with size N, we calculated the pairwise correlations among cell types of tissues and identified maximally and minimally correlated cell types in the 3-months age group. Then, we calculated age-related changes of those correlations using Spearman correlation to construct null distribution.

Analysis of within tissue convergence of cell types
Analogous to inter-tissue convergence analysis, we also studied intra-tissue convergence of cell types in scRNA seq data by calculating CoV among cell types within a tissue for each individual in 3m, 18m and 24m separately. We filtered the data to obtain cell types present in at least 2 individual mice in every time point for each tissue which yielded 4, 7, 20 and 6 cell types in brain, liver, lung and muscle, respectively. We then tested the mean CoV (or CoV per gene) change with age using Spearman correlation test.

Ethics Statement
Post-mortem samples were obtained from 16 C57BL/6J mice aged between 2 days and 904 days, housed in groups and fed under standard conditions in the MPI-EVA Animal Facility (Leipzig, Germany). The mice were sacrificed for reasons independent of this study, their tissues were 25 harvested and frozen immediately, and stored at -80°C.

Competing Interests
The authors report no competing interests.

Data Availability
Raw and processed RNA-seq data have been deposited in GEO with the accession number GSE167665. All summary statistics and analysis output are also provided as supplementary tables.

Code Availability
All the code used to perform analyses is available in GitHub: https://github.com/hmtzg/geneexp_mouse   (Figure 1-figure supplement 4) (Figure 1-source data 1).   Figure 1-source data 1.        Figure 2-source data 1).     The number of CoV changes with age (without a significance cutoff) during development and ageing.

GORA of age-related genes in tissues
Tissue-specific age-related gene expression changes and functional enrichment test results performed with gene over-representation analysis (GORA) using 'topGO' package.

Supplementary File 2. GORA of shared age-related genes among tissues
Functional enrichment for shared genes across tissues. The same GORA that was performed for Supplementary File 1, was used to test the enrichment of shared up/down-regulated genes in development among the background genes which are chosen as the all significant age-related genes across tissues in development. We did not apply the test for the ageing period as there were no shared ageing-related expression changes ( Figure   S5).

Supplementary File 3. GORA of reversal patterns
Functional enrichment for gene expression reversals.GORA analysis was performed with the same criteria as explained above. Up-Down reversal genes were tested against Up-Up genes and Down-Up reversal genes were tested against Down-Down genes in each tissue.

Supplementary File 5. GORA of DiCo gene clusters determined with expression levels
Functional enrichment of DiCo genes clustered with kmeans algorithm according to their expression levels. Gora analysis was performed using gene sets in each cluster (Figure 2-figure supplement 4) which are tested among all DiCo genes.