Characterizing expression changes in noncoding RNAs during aging and heterochronic parabiosis across mouse tissues

Mapping of ncRNA expression across mouse organs

We sequenced 771 tissue samples of the TMS cohort to map molecular shifts across the whole organism during healthy aging (Fig. 1a). The protocol enriches for small ncRNA, especially mature miRNAs. Even though full-length reads cover only small ncRNAs (miRNAs or piwi-interacting RNAs (piRNAs)) completely, the protocol generates measurable fragments of longer ncRNAs. This sequencing strategy extends the existing mRNA TMS dataset⁵ with miRNA, piRNA, long ncRNA (lncRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), transfer RNA (tRNA), ribosomal RNA (rRNA) and small Cajal body-specific RNA (scaRNA). The tissue sample collection includes 16 solid tissues of C57BL6/JN mice (bone (femurs and tibiae), brain (hemibrain), brown adipose tissue (BAT, interscapular depot), gonadal adipose tissue (GAT, inguinal depot), heart, kidney, limb muscle (tibialis anterior), liver, lung, bone marrow, mesenteric adipose tissue (MAT), pancreas, skin, small intestine (duodenum), spleen and subcutaneous adipose tissue (SCAT, posterior depot)). The selected time course covers the mouse lifespan from a developmental age of 1 month up to 27 months (males: aged 1, 3, 6, 9, 12, 15, 18, 21, 24 and 27 months; females: aged 1, 3, 6, 9, 12, 15, 18 and 21 months). With up to six mice per timepoint, the study covers a maximum of 960 samples (16 organs × 10 timepoints × 6 replicates). As not all mice survived to the later timepoints and we further excluded 26 low-quality RNA samples, we finally included 771 high-quality samples in the study (Supplementary Table 1).

a, Study overview—data of the aging (TMS) cohort, consisting of mouse samples collected from 16 different tissues at ten different timepoints throughout the lifespan with maximal six replicates per timepoint varying due to sample and sequencing quality. A total of 771 samples were sequenced, and the reads were annotated to the 87,590 different RNA reference sequences from eight RNA classes. Of the RNAs in the databases listed on the left, 58,422 different ncRNAs were annotated in the raw reads and we found 7,883 noncoding features as abundant expressed in our TMS aging cohort. Created with BioRender. b, t-SNE visualization of all samples of the TMS cohort over all detected noncoding RNAs, colored by tissue of origin. c, Percentage of counts per RNA class, calculated on total counts per tissue after local filtering for all tissues in the TMS cohort, color coded by RNA-class color legend as indicated in a. d, Variation of mean count distribution per RNA class over the lifespan of the mouse in the brain, BAT, marrow and liver; calculated count percentages per sample after local filtering. Created with BioRender. e, Density plot of Spearman rank correlation of all expressed noncoding RNA with age in each individual tissue grouped by RNA classes, density scaled individually for every RNA class.

We mapped resulting sequencing reads against 87,590 ncRNA sequences (Fig. 1a, left column) derived from established reference databases (miRNAs, miRBase 22, tRNAs: GtRNAdb 18.1, piRNA: RNACentral 15, all other ncRNAs: Ensembl 100). Altogether, we detected reads mapping to 58,422 different ncRNAs (Fig. 1a, middle column), with miRNAs being the most abundant class. An average of 36.2% of reads across tissues mapped to miRNAs (Extended Data Fig. 1a). The distribution of reads to RNA classes, however, varied substantially between tissues (P < 0.05, Kruskal–Wallis test; Extended Data Fig. 1b). We thus asked whether the variation in read distribution is related to the length of representatives. We generated aligned sequence profiles to quantify the length of sequences covered with our reads. First, we explored the percentage of sequence length covered versus the sequence reference length. Even for very long sequences exceeding 10,000 bases, we partially recovered large fractions or even the complete sequence (Extended Data Fig. 1c). We also computed the maximal assembly for each RNA, that is, the longest contiguous read mapping, and compared it with the sequence reference length (Extended Data Fig. 1d). Even though the fraction of the sequence covered by the maximal assembly decreased for larger RNAs (Spearman’s rho = −0.43), we verified that throughout all RNA classes, we were still able to reproduce for a subset up to 100% of the full-length reference (Supplementary Table 2). Therefore, we decided to include the full dataset as a reference for future studies while implementing restrictive filtering steps to increase the reliability of our data. Especially the somatic piRNAs exceeded the expected counts, likely driven by artifacts in piRNA annotation¹³ and calling for additional quality control filters. We first retained piRNAs encoded in prepachytene piRNA genomic clusters^14,15 to minimize the number of false positive hits. Next, we removed low abundant features across all ncRNA classes, keeping those with at least 1 read mapped per million (rpmm) in at least one sample, resulting in the abundant dataset (Fig. 1a, right column). Applying this stringent filtering, the number of piRNAs in our dataset decreased from 43,799 detected down to 43 abundant, likely removing most falsely annotated features¹³.

Clustering by ncRNA expression using t-distributed stochastic neighbor embedding (t-SNE), samples split into tissue-specific groups (Fig. 1b). One cluster contained skin, GAT and SCAT samples, which likely can be explained by their biological and functional close relationship of containing similar cell types. To check whether relevant biological factors outweigh technical ones, we performed a principal variance component analysis. The highest proportion of variance in the data was explained by tissue identity (Extended Data Fig. 2a). Annotating the t-SNE plot by animal sex revealed a uniform spread, excluding it as a major driver of the observed variance (Extended Data Fig. 2b).

Following our main objective to identify organ-specific aging trajectories, we added a tissue-specific, that is, local filtering to check whether ncRNA expression changed not only between tissues but also with age (cf. Methods). On the locally-filtered data, we calculated read count percentages for all RNA classes. As for the detected reads, we observed tissue-specific distributions (Fig. 1c). Analyzing those over time, we identified two clusters of tissues (Extended Data Fig. 2c). One exhibited a stable count distribution (mean variance < 4.5%) and the other showed high variance within the count distribution (mean variance > 4.5%). Specifically, 3 of 16 tissues showed high variance (brain, BAT and limb muscle), while most tissues (13–16) were characterized by a stable read distribution (including marrow and liver) (Fig. 1d and Extended Data Fig. 3). In the brain, the share of snRNA reads decreased from 77.9% to 10.0%, while the share of miRNAs increased from 9.1% to 28.5%. In BAT, the miRNA share grew steadily from 4.1% to 26.4% and the rRNA share dropped from 62.7% to 27.8%.

The observed variations of the RNA classes prompted us to assess the expression changes during aging for the individual ncRNAs. Therefore, we determined the Spearman rank correlation of age with the expression of every ncRNA in each tissue separately. We identified 31 tRNA fragments that were substantially differentially expressed between 3 and 21 months (two-sided t-test, P adjust < 0.05). Eight tRNA fragments showed increased expression (in brain and lung) and 23 showed decreased expression with age (in bone, limb muscle, skin and GAT). tRNA-related metabolism, transcription, modification and derivatives have vital roles in aging and longevity of organisms, as tRNA expression decreases with age¹⁶. We further observed that miRNAs displayed the strongest correlations with age over all tissues (exceeding the interval of −0.5 to 0.5; Fig. 1e). Given that miRNAs were captured in full-length by our sequencing platform, their high abundance across tissues and the fact that they exhibited the largest effect size, we further focused on miRNAs for downstream analysis.

MiRNA lifespan trajectories are largely tissue specific

For the intersection of miRNAs expressed in all tissues, we observed more markers being correlated positively than negatively with age (Fig. 2a and Supplementary Table 3). In contrast, large sets of miRNAs were correlated with age in a specific tissue. For example, six miRNAs were negatively correlated exclusively in limb muscle and 37 were positively correlated only in BAT. One of these miRNAs, miR-107, regulates insulin sensitivity and is postulated as a target for the treatment of type 2 diabetes and obesity¹⁷. Its increase in aging could be connected to the fact that age is a risk factor for diabetes¹.

a, Heatmap of Spearman rank correlation values of the intersection of miRNAs expressed in all tissues, color coded for positively correlated in blue (r > 0.5), negatively correlated in red (r < −0.5) and not correlated in white (−0.5 < r < 0.5). b, Heatmap of miRNAs (anti-) correlated with age in at least two tissues, colored by number of tissues (anti-) correlated and divided into miRNA positively and negatively correlated with age. c, Heatmap for the count of deregulated miRNAs in each tissue at each subsequent timepoint. Deregulated miRNAs are determined by calculating the foldchange of all later timepoints versus 3 months of age and miRNA with foldchanges <2/3 or >3/2 are considered deregulated. d, Volcano plot of all miRNAs expressed in brain, log₂(FC) versus −log₁₀(P values) (two-sided t-test) calculated between mice aged 3 months and all later timepoints with comparisons for 12 months (light red) and 18 months (light blue) highlighted. e, Whole organism miRNA trajectory clustering—z-scored trajectories of each expressed miRNA in each tissue over the entire lifespan of the mice were calculated. These trajectories were grouped into 20 clusters. Three clusters are displayed as examples, showing tissue-specific miRNA signatures. Cluster 2 is composed of 50.1% miRNAs originating from skin, cluster 9 of 61.2% miRNAs from brain and cluster 13 of 70.6% miRNAs from BAT. f, Cluster 20 of the whole organism miRNA trajectory clustering—the cluster contains two global aging miRNA, miR-29a-3p from ten different tissues and miR-29c-3p from eight different tissues.

Certain miRNAs were linearly correlated with age in more than one tissue (Fig. 2b). MiR-29a-3p was positively correlated in eight tissues, and miR-300-3p, miR-487b-3p and miR-541-5p were negatively correlated in five tissues each. Based on these observations, we separated the miRNAs into the following three classes: nonaging-related, local aging and global aging miRNAs. Local aging miRNAs were defined as correlated with age in at least one tissue exceeding the interval of −0.5 to 0.5. We accordingly defined miRNAs correlated with age in more than five different tissues as globally aging. Following these definitions, we identified the three mentioned negatively correlated miRNAs together with five positively correlated miRNAs (miR-29a-3p, miR-29c-3p, miR-155-5p, miR-184-3p and miR-1895) as globally aging.

We then examined whether nonlinear age-related expression changes occur as well. Using the 3m timepoint as a baseline, we calculated foldchanges (FC) for all later timepoints and respective P values. Based on this analysis, we determined the number of deregulated miRNAs (Fig. 2c). Most were deregulated in BAT, driven by the large fraction of positively correlated local aging miRNAs. Investigating the brain, we found a peak at the ages 12 and 18 months with a count of 412 and 427 deregulated miRNAs, respectively. Most of all substantially deregulated miRNAs (77.6%) in brain showed the strongest effect at 12 and 18 months (Fig. 2d, Extended Data Fig. 4a,b and Supplementary Table 4). The higher count of deregulated miRNA at certain timepoints matched our expectation, as we hypothesized that miRNAs were responsible for the regulation of the previously reported transcriptome changes¹⁸. We further confirmed that those effects were not driven by lowly expressed features—we projected the mean expression against the FC for all tissues and all ncRNAs per timepoint (Extended Data Fig. 5). In line with our assumption, substantial FCs could be observed across all expression scales.

To identify common patterns within the nonlinear changes over time, we calculated z scores for all miRNAs being expressed in every tissue. Each miRNA in every single tissue was displayed as an aging trajectory and clustered across all organs. Ten of the 20 clusters obtained were composed mainly of one tissue; thus, we propose the existence of organ-specific miRNA time course signatures (Fig. 2e). Half of the miRNAs in cluster 2, with a peak at 3 months and a late increase again at 24 months, originated from the skin. Cluster 9, which showed a peak at 12 and 18 months, was composed of 61.2% brain miRNAs. The expression of miRNAs in cluster 13 increased continuously from the age of 6 months on, and this trajectory was specific for BAT (70.6%). In summary, we determined 10 of the 20 clusters to be tissue-specific, with at least 30% of miRNA originating from a single tissue (Extended Data Fig. 6).

The global aging miRNAs marked an exception to this tissue-specific clustering. Trajectories from more than five different organs for seven global aging miRNAs clustered together. For instance, we found the trajectories of miR-29a-3p and miR-29c-3p from ten and eight different tissues in cluster 20, respectively (Fig. 2f). The expression of miRNAs within this cluster increased continuously with age. This consistent signature could be indicative of the regulation of key pathways across all organs upon aging. Thus, we investigated the relationship between miRNA and mRNA expression closer.

Transcriptome changes mirrored by global aging miRNAs

The previous analyses suggested five miRNAs as cross-organ aging markers increasing with age (Fig. 2b). Following the biological mechanism, we expected repression of target genes with aging. We chose to identify potential new targets in an unbiased manner by correlating miRNA with mRNA expression levels from the TMS dataset⁵. In the first step, we defined targets by exhibiting a significant inverse correlation (r < −0.4, P < 0.05). To support the validity of our approach, we checked the share of predicted miRNA–mRNA interactions with conserved binding sites for the miRNAs. As a control, we compared this number against the share of conserved binding sites in the miRNA–mRNA interactions predicted via positive correlation. For 7.3% (9 of 122) of the miRNA–mRNA interactions identified via inverse correlation, we found at least one conserved miRNA binding site, as compared to the 2.1% (120 of 54,992) miRNA–mRNA interactions in the control set (Fisher’s exact test, P = 0.0018). Because a gene can contain multiple binding sites across multiple 3′ UTRs and different site types exhibit different strengths, we repeated the analysis for each type of binding site. The amount of conserved 8mer binding sites is 6.3 times higher as compared to the control (4.91% inv. correlation, 0.78% control; P = 0.0006), for conserved 7mer-8m binding sites 3.8 times higher (4.92% inv. correlation, 1.30% control; P = 0.0062) and for conserved 7mer-1a binding sites 9.0 times higher (2.45% inv. correlation, 0.27% control; P = 0.0064).

The filtered target gene sets showed distinct overlaps (cf. Methods; Fig. 3a and Supplementary Table 5). Three of the six targets are shared among all miRNAs, Eln, Col1a1 and Col3a1, which have a role in protein digestion and absorption and encode extracellular matrix (ECM) proteins. These are already validated targets for miR-29b-1/miR-29a (ref. ¹⁹). Overall, enriched processes for all targets were dominated by ECM-associated processes, such as ECM organization, collagen fibril organization and ECM-receptor interaction (Fig. 3b)²⁰. Senescent cells are known to exhibit altered expression and organization of ECM and the ‘senescence-associated secretory phenotype’^1,21. Our data suggest that these effects could be regulated by global aging miRNAs. Another part of the network composed of mainly Y-chromosome-coded proteins contained proteins related to ‘ubiquitin-proteasome dependent proteolysis’ (Usp9y), histone modification introducing proteins (Kdm5d) and probable transcriptional activators (Zfy1, Zfy2). Hence, other layers of regulation mechanisms are targeted. The ‘AGE-RAGE signaling pathway in diabetic complications’ and ‘dysregulated miRNA targeting in insulin/PI3K-AKT signaling’ were enriched, supporting our suggestion of the importance of miRNA regulation in nutrient sensing.

a, Venn diagram of predicted target transcripts of the five global aging miRNAs positively correlated with age in most tissues. Targets are identified via inverse correlation of expression values (Spearman’s rho < −0.4, Spearman’s statistics P < 0.05, two-sided); only miRNA–mRNA target predictions were selected that are correlated in at least two tissues for one of the five miRNAs. b, STRING network for all connected proteins encoded by target transcripts of the global aging miRNAs positively correlated with age (a); nodes are color coded for pathways in red for ‘ECM’, purple for ‘secreted’ and green for ‘dysregulated miRNA targeting in insulin PI3K-AKT signaling’. c, Venn diagram of predicted target transcripts of the three global aging miRNAs negatively correlated with age in most tissues. d, STRING network for all connected proteins encoded by target transcripts of the negatively with age correlated global aging miRNA negatively correlated with age (c); nodes are color coded for pathways in red for ‘immune receptor activity’, in purple for ‘cytokine activity’, in yellow for ‘hematopoietic cell lineage’, in pink for ‘adaptive immunity’ and in green for ‘immunoglobulin’. e, Expression of miR-29c-3p in reads per mapped million in the liver (r = 0.69), kidney (r = 0.56), lung (r = 0.51) and BAT (r = 0.48) over the mouse lifespan (mean per timepoint ± s.d.), created with BioRender.

Consistent with the enriched pathways for the targets of cross-organ aging miRNAs were the enriched pathways for the targets of the local aging miRNAs. The ‘PI3K-AKT signaling pathway’, ‘protein digestion and adsorption’, ‘metabolic pathways’, ‘adipocytokine signaling pathway’ and ‘insulin resistance’ were found among the top 20 locally enriched pathways in targeted mRNAs (Extended Data Fig. 6b).

Through a reduction of miRNA expression during aging, the repression of gene expression is potentially reduced or even lost. Gene targets for global aging miRNAs reducing repression with age were identified via correlation (Fig. 2b) (r < −0.4, P < 0.05). The expression of three miRNAs, miR-300-3p, miR-487b-3p and miR-541-3p, decreased during the lifespan in five tissues. The overlap of their potential targets was high, with 138 of 327 predicted interactions (Fig. 3c and Supplementary Table 6). The identified targets exhibited a functional enrichment for pathways related to immune system processes, such as ‘cytokine–cytokine receptor interaction’, ‘Th1 and Th2 cell differentiation’, ‘Th17 cell differentiation’, ‘chemokine signaling pathway’ and ‘NF-kappa B signaling pathway’. The network was particularly dense in its center, with targets related to ‘adaptive immunity’, ‘immunoglobulin’, ‘hematopoietic lineage’, ‘immune receptor activity’ and ‘cytokine activity’ (Fig. 3d). We also determined the locally enriched pathways for all miRNAs in every individual tissue whose expression decreased upon aging via inverse correlation with mRNA targets. These were similarly dominated by immune-related processes (Extended Data Fig. 6c). As immune senescence and inflammation are hallmarks of aging¹, it is crucial to further investigate these potentially age-sensitive regulation mechanisms.

We chose the global aging miR-29c-3p as an example for further investigation. In liver and kidney, expression increased monotonically over the lifespan as well as in BAT but at a lower baseline expression (Fig. 3e). In the lung, the steep increase during early adulthood ends at approximately 12 months of age. A general trend of miR-29c-3p expression increase was present in all tissues, but expression levels and the course of increase showed tissue-specific patterns (Extended Data Fig. 7).

miR-29c-3p exhibits an organ-specific rejuvenation response

Expansive beneficial effects on cognition, muscle strength and bone repair have been observed for heterochronic parabiosis via a shared common circulation, or systemic infusions of young blood²². We sequenced tissue samples from a parabiosis intervention cohort to determine whether the young blood in aged individuals influences small ncRNA expression. The cohort was composed of 176 samples from six different organs of isochronic young (IY) and aged (IA), and heterochronic young (HY) and aged mice (HA) (Supplementary Table 1). Rejuvenation, the reversion of aging aspects, is the desired outcome of the intervention. However, it is accompanied by accelerated aging, the negative effect of the young sharing their blood with the old. In our study, the rejuvenation effect was measured by comparing the expression levels in IA mice with those detected in HA mice. In turn, the accelerated aging effect was defined by the difference between IY and HY mice (Extended Data Fig. 8a). Healthy aging was defined as the comparison of mice from the TMS cohort aged 3 and 21 months (AGE), closely matching the age distribution of the parabiosis cohort at takedown. Clustering the samples using t-SNE revealed tissue identity as major driver of variance across the experimental groups (Extended Data Fig. 8b,c).

We assigned deregulated miRNAs to the following groups: either (1) uniquely deregulated in rejuvenation (REJ unique) or in accelerated aging (ACC unique), or (2) deregulated in physiological aging as well as rejuvenation (REJ up and AGE down, REJ down and AGE up) or accelerated aging (AGE and ACC up/down). We found 233 uniquely deregulated miRNAs in rejuvenation and 43 in accelerated aging (Fig. 4a). Intriguingly, 17 age-related miRNAs were deregulated in the opposite direction in REJ. No miRNAs were deregulated in AGE and in the same direction in ACC, but the uniquely rejuvenated miRNAs were enriched in certain pathways in MAT (‘insulin resistance’, ‘adipocytokine pathway’, ‘type 2 diabetes mellitus’), which again have a role in nutrient sensing.

a, Heatmap of the absolute numbers of deregulated miRNAs in parabiosis for the six comparisons in all analyzed tissues—uniquely deregulated in REJ, uniquely deregulated in ACC, in REJ up-regulated and healthy aging (AGE) downregulated, and vice versa, in AGE and ACC up-regulated or downregulated, all values >100 set to 100. b, Boxplot of expression z-scores for miR-29c-3p per timepoint in the liver for healthy aging from the TMS cohort and z-scores of miR-29c-3p in the liver from the parabiosis cohort for the following four different groups: IY, HY, HA and IA; arrows on the right indicate the effect size of aging in the parabiosis cohort, rejuvenation and accelerated aging. Boxes span the first to the third quartile with the median value represented as line inside the box. The whiskers show the maximum and minimum values or values up to 1.5 times the interquartile range above and below the first or third quartile if outliers are present. c, Venn diagram of local aging miRNAs (exclusively positively correlated with age in one tissue; green), circulating miRNAs (all miRNAs expressed in EVs and serum, red) and global aging miRNAs (positively correlated with age in more than five tissues, blue). d, Summary of validated targets of miR-29 family members, shown in their cellular location. These targets are key proteins in processes related to the ECM; newly validated targets are shown in red. Created with BioRender.

For three global aging miRNAs, we discovered that changes in expression observed during healthy aging can be partially reversed in response to parabiosis. For miR-29c-3p, we measured a strong rejuvenation effect in the liver, four times higher than the effect of accelerated aging (Fig. 4b). The other two global aging miRNAs miR-184-3p in the liver and miR-300-5p in GAT showed similar trends of reversed expression but with a lower magnitude (Extended Data Fig. 9a,b). Considering the pronounced globally aging versus local rejuvenation profile of miR-29c-3p, we chose to explore systemic effectors and mediators of these signals.

Expression of circulating mir-29 family increases with aging

MiRNAs can circulate in the plasma and EVs between organs. We thus assessed the abundance of miR-29c-3p in both plasma and the vesicle-bound fraction using an independent cohort²³. Analyzing the expression at five timepoints across the lifespan from 2 to 18 months allowed us to correlate and compare the abundance of the miRNA in plasma and EVs. We observed an increase of miR-29c-3p expression correlated with age for both fractions (r = 0.56 (plasma) and 0.65 (EVs)). The share of positive global aging miRNAs detected as circulating was higher than the share of local aging miRNAs (38.3%) (Fig. 4c). This supports our hypothesis that miRNAs traveling via the shared circulatory system could have a role in the positive effect of parabiosis.

MiR-29c-3p could regulate gene expression in pathways resulting in health improvements by entering the tissue via the blood in vesicles. Recently, miRNAs with certain sequences were shown to be more likely secreted in small EVs, and their capability to inhibit target genes in recipient cells is enhanced⁸. One so-called EXOmotif is CNGGNC, which is very similar to a sequence found in the mature mmu-miR-29c-3p CUGGUG. We performed luciferase assay experiments to validate our predictions for mir-29-family members on target genes related to aging. Lox and Adamts17 were validated as high confidence targets, and Vash1 was validated as a low confidence target (Extended Data Fig. 9c,d). Previously known targets from the literature (Eln, Col1a1, Col1a2, Col3a1 and Adam12), as well as Lox and Adamts17, are components in ECM processes (Fig. 4d), supporting our hypothesis that mir-29-family members have a crucial role in organismal aging due to their repressive regulatory function on these targets.