Research Paper Volume 13, Issue 14 pp 18150—18190
Integration of segmented regression analysis with weighted gene correlation network analysis identifies genes whose expression is remodeled throughout physiological aging in mouse tissues
- 1 Institute of Biomedicine – iBiMED, Department of Medical Sciences, University of Aveiro, Aveiro 3810-193, Portugal
- 2 Life and Health Sciences Research Institute (ICVS), School of Medicine, University of Minho, Braga 4710-057, Portugal
- 3 ICVS/3B’s–PT Government Associate Laboratory, Braga/Guimarães, Portugal
Received: February 19, 2021 Accepted: July 21, 2021 Published: July 29, 2021https://doi.org/10.18632/aging.203379
How to Cite
Copyright: © 2021 Ferreira et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Gene expression alterations occurring with aging have been described for a multitude of species, organs, and cell types. However, most of the underlying studies rely on static comparisons of mean gene expression levels between age groups and do not account for the dynamics of gene expression throughout the lifespan. These studies also tend to disregard the pairwise relationships between gene expression profiles, which may underlie commonly altered pathways and regulatory mechanisms with age. To overcome these limitations, we have combined segmented regression analysis with weighted gene correlation network analysis (WGCNA) to identify high-confidence signatures of aging in the brain, heart, liver, skeletal muscle, and pancreas of C57BL/6 mice in a publicly available RNA-Seq dataset (GSE132040). Functional enrichment analysis of the overlap of genes identified in both approaches showed that immune- and inflammation-related responses are prominently altered in the brain and the liver, while in the heart and the muscle, aging affects amino and fatty acid metabolism, and tissue regeneration, respectively, which reflects an age-related global loss of tissue function. We also explored sexual dimorphism in the aging mouse transcriptome and found the liver and the muscle to have the most pronounced gender differences in gene expression throughout the lifespan, particularly in proteostasis-related pathways. While the data showed little overlap among the age-dysregulated genes between tissues, aging triggered common biological processes in distinct tissues, which we highlight as important features of murine tissue physiological aging.
Gene expression alterations occurring throughout the lifespan have been described for a multitude of species, organs, and cell types [1–10]. The most commonly reported age-related dysregulations involve the immune system [9, 11–13] where inflammatory response genes are upregulated even in the absence of pathogen infection [5, 6, 9, 11, 14–19]. Energy metabolism, redox homeostasis, and mitochondrial function alterations are also frequently observed in age-related studies [6, 9, 11, 15–18, 20], particularly the downregulation of genes encoding mitochondrial ribosomal proteins and components of the electron transport chain [5, 11, 14–16, 18], protein synthesis machinery [5, 11, 17], developmental and cell differentiation genes [9, 11, 19], and extracellular matrix components [6, 14–16]. Up-regulated genes are associated with the stress response and DNA repair [5, 6, 9, 11, 14, 16–18], RNA processing [11, 12, 17] and cell cycle arrest [5, 16, 19]. Despite this, the existence of specific genetic signatures of aging continue to be a matter of debate as gene regulation is mostly tissue- [5–7, 15, 20–24] and cell-specific [25, 26], but also because there is focus on comparisons between young and old individuals without much consideration of the dynamics of gene expression throughout the lifespan. Nonetheless, recent studies in humans and animal models shed light on these dynamics. As an example, a marked shift in mRNA and microRNA expression has been reported to occur at around age 20 in the human prefrontal cortex [27, 28]. Less striking alterations were reported in the same brain region between 30-60 years [28, 29], entailing genes related to the synapse, fatty acid metabolism, purine nucleotide binding, ubiquitin proteolysis, channel activity, translation, DNA damage response, transcriptional activation, and neuronal function . Late middle-age and early old-age shifts have also been described in human peripheral blood leukocytes for genes involved in cancer, hematological and immunological diseases, cell-mediated immune response and signaling pathways , and in the human brain and muscle for both coding and non-coding RNAs pertaining to longevity pathways .
Similar findings have been reported in a study across 11 rat organs where the most frequent changes in gene expression occurred at around 6 and 21 months , proposed to be equivalent to middle-age in humans . Another study across 17 mouse tissues, whose dataset we re-analyzed in this work, identified shift points of gene expression trajectories at around 6 months for extracellular matrix genes, 10 months for mitochondrial genes, 12 months for genes encoding heat shock proteins, and at around 15 months for immune response genes . Other studies also suggest tissue-specific turning points in gene expression profiles [6, 20, 35]. For example, immune response gene expression was found to change in the mouse kidney between 13 and 20 months, in line with the previously described organismal trend, whereas in the spleen and lung this shift occurs later in life, at around 26 months .
Interestingly, various works show prominent sex-differences in gene expression and regulation in mammals, affecting processes such as hormone secretion, immune response, extracellular matrix organization, oxidoreductase activity, lipid metabolism, nucleotide metabolism, cytoplasmic and mitochondrial translation, RNA helicase activity, ribosomal RNA processing, synaptic plasticity, and neurotransmitter transport [36–38, recently reviewd in 39]. However, evidence regarding gender-biased gene expression across the lifespan is much scarcer. A recent study in mice found sex-specific differences in the aging of cells of renin lineage in the kidney, with aged females exhibiting up-regulation of genes involved in angiogenesis, apoptosis, epithelial to mesenchymal transition, and TGFβ signaling, whereas in aged males these genes were down-regulated . Another study showed age-related induction of genes involved in the apoptosis of microglia in the retinas of old female mice, while in old males the expression of these genes was not affected by aging . Sex differences across the lifespan have also been described for the mouse hippocampus, with aged females exhibiting activation of genes involved in inflammatory processes, when compared with aged males .
These studies have mainly focused on differential gene expression to characterize the aging transcriptome, using methods that rely on static comparisons of mean expression values between consecutive age groups, or relative to a reference time point, thus ignoring the influence of gene expression levels in earlier time points in those of later ages [comprehensively reviewed in 43] and making it difficult to perform gene prioritization. These methods also disregard the pairwise relationships between gene expression profiles, which may underlie commonly altered pathways and regulatory mechanisms with age. To clarify these issues, we have combined segmented regression analysis with weighted gene correlation network analysis (WGCNA) to re-analyze a publicly available mouse RNA sequencing (RNA-Seq) dataset (GSE132040) [6, 44].
The GSE132040 dataset was made publicly available by the Tabula Muris Consortium  and consists of transcriptomic data from 17 male and female mouse tissues sampled across 10 different time points, making it a very useful resource for the study of mammalian aging. In the article resulting from the analysis of this dataset , Schaum and colleagues performed differential gene expression, both relative to a reference (1, 3 or 6 months) or between consecutive time points and, acknowledging the limitations of pairwise comparisons, also grouped genes into different clusters based on Euclidean distance and fitted local regression (LOESS) models for each gene in order to address the expression dynamics across the lifespan. We believe this approach does not fully account for the limitations we presented above, and since our goal is to identify gene expression signatures and shifts in expression throughout aging, we decided to re-analyze this data using the R packages Trendy  and WGCNA .
Trendy is an R package that fits a set of segmented (piece-wise) regression models to high-throughput, ordered expression data, allowing to identify dynamic gene expression patterns over a time course and to determine the direction and point in time when changes in expression occur . It is a recently developed methodology that has been so far applied in stem cell differentiation  and in hepatic disease progression  studies. Interestingly, the latter research also integrated Trendy with WGCNA to identify high-confidence genes involved in human autoimmune hepatitis (AIH) . Conversely, WGCNA is a widely implemented systems biology approach not only used for constructing co-expression networks based on the pairwise correlations between genes, but also for detecting significant associations between these networks’ modules and a given trait of interest, and to identify important players in these associations . It has been broadly applied in the context of aging research to identify novel age-related gene expression signatures, as is the case of CR2, VPREB3, MS4A1 and CCR6, involved in B cell activation and receptor signaling pathways in human peripheral blood cells , PPP3CB, CAMSAP1, ACTR3B, and GNG3, involved in synaptic vesicle cycle, cGMP-PKG, and dopaminergic synapse signaling pathways in human pre-frontal cortex , and PGLS, a gene potentially involved in synaptic loss in the aging brain of rhesus macaques , among other examples.
For all these reasons, we subjected the Tabula Muris Consortium’s GSE132040 dataset to a new analysis pipeline comprising the Trendy and WGCNA approaches in order to: 1) identify expression signatures of genes significantly correlated with age (and with sex) in different mouse tissues, 2) establish their trajectories from mature adulthood to old age, 3) identify the time point of the shift in the gene expression profile, and 4) evaluate the biological processes (BPs) associated with the gene dysregulations. Different onsets of gene dysregulation were identified, demonstrating the asynchronous impairment of gene expression with age described by Schaum et al. . Gene Ontology (GO) biological processes’ over-representation analysis revealed that inflammation and immunity-related responses are prominently altered in the brain and the liver, while in the heart and the muscle aging affects fatty acid metabolism, and tissue regeneration, respectively, which may reflect the global loss of organ function, confirming previous reports. Additionally, functional enrichment analysis showed sex-differences in the expression throughout the lifespan of genes involved in energy metabolism and proteostasis-related pathways, with only males exhibiting age-related dysregulation of these processes. Despite little overlap in age-dysregulated genes between tissues, a comparison at the level of dysregulated processes revealed inter-tissue commonalities, such as alterations in immune response, tissue regeneration, energy metabolism, glucocorticoid signaling and response to amino acid stimulus. We propose that the genes involved in these processes are important players of murine physiological aging.
Most dynamic changes in gene expression across the lifespan are observed within mature adulthood and middle life
In order to obtain a global characterization of the mouse aging transcriptome and, in particular, to understand which of the main known sources of genetic variation (tissue, age, and sex) is responsible for the highest percentage of sample segregation, we performed a principal component analysis (PCA) of brain, heart, liver, muscle and pancreas samples, ranging from 3- to 27-month -old male and female mice, based on variance stabilizing transformation (VST)-normalized read counts. Principal components were calculated based on the 500 most variable genes as they are expected to capture the greatest variability between samples (see Methods - Data pre-processing and normalization. We observed that, based on VST-normalized gene expression values, samples tend to cluster by tissue, which is in line with what Schaum et al.  found, as well as with previous observations [5, 7, 15, 20–24]. For this reason, subsequent analyses were performed on each tissue separately (Figure 1A). Additionally, we focused only on samples from the brain, heart, muscle, liver, and pancreas, in order to provide an in-depth discussion of the findings and also for future integration with in-house proteomics data. When considering gene expression variation in each tissue independently, we observed that sex, rather than age, is responsible for most of the between-sample variability and adjusted for its effect on gene expression by adding it as a co-variable in the regression model (Supplementary Figure 1; see Methods – Trendy segmented regression analysis.
Figure 1. Whole-transcriptome characterization of different mouse tissues throughout the lifespan by segmented regression analysis. (A) PCA of all tissues performed on VST-normalized read counts of the 500 most variable genes and colored by all known effects highlights type of tissue as the main contributor to sample segregation. (B) Percentage of top dynamic (Trendy) genes across the lifespan for the mouse brain, heart, liver, muscle, and pancreas, highlighting the liver and the pancreas the ones with higher and lower dysregulation, respectively. The displayed values correspond to the percentages of the total number of Trendy genes found in each tissue relative to the total number of expressed genes per tissue (brain: 34164, heart: 28073, liver: 20157, muscle: 18978, pancreas: 18414). Top dynamic genes were selected based on tissue-specific adjusted R2 thresholds (brain: > 0.2; heart: > 0.2; liver: > 0.1; muscle: > 0.3; pancreas: > 0.3) and p-values < 0.1 in at least one segment; Supplemental File SF1). (C) Histograms of the distribution of breakpoints in gene expression of the top dynamic genes per tissue. Each bar depicts the sum of all breakpoints of all Trendy genes at that given time point. Monotonic behaviors (i.e., no breakpoints) are not included in the histograms. Dotted, vertical lines indicate the median breakpoint distribution for each tissue. Median breakpoint distribution in the muscle was significantly different from all the other tissues (Kruskal-Wallis Test followed by Dunn’s Test; Supplementary Table 1). (D) Percentage of top dynamic genes exhibiting 0 to 7 breakpoints (maximum number of breakpoints allowed in the Trendy regression model). In all tissues, except from the Pancreas, most Trendy genes exhibit monotonic expression patterns (continuously up or down). In the Pancreas, the majority of genes display one breakpoint. (E) Biotype distribution of the top dynamic genes per tissue. Biotype nomenclature based on Ensembl annotation. In all tissues, protein coding genes were significantly over-enriched relative to the reference genome’s annotation, whereas lncRNAs were significantly under-enriched (Fisher’s Exact Test; Supplementary Table 2).
As a first approach to establish age-regulated genes, we carried out segmented regression analysis on each tissue’s expressed genes (brain: 34164, heart: 28073, liver: 20157, muscle: 18978, pancreas: 18414; post-filtering, normalized expression; see Methods – Data pre-processing and normalization), using the R package Trendy (v. 1.8.2) . In brief, in a segmented regression model, each gene’s expression is regarded as a linear piece-wise function over time, with each segment being separated by a breakpoint (i.e., a point when gene expression changes). The Trendy model fits to each gene multiple segmented regression models with varying number of breakpoints, and then selects the one who explains the best the dynamics of that given gene’s expression across the time-course (see Methods – Trendy segmented regression analysis). Then, we selected the top dynamic (Trendy) genes based on the goodness of fit of the chosen model to the corresponding gene expression pattern, and on the significance of the segments (see Methods – Trendy segmented regression analysis; Supplementary Table 10).
Next, to obtain a global view of the amplitude of gene expression dysregulation of each tissue with aging, we calculated the percentage of Trendy genes across the lifespan, considering each tissue’s total number of top dynamic genes relative to each tissue’s total expressed genes. We identified a total of 747, 1799, 3690, 1048, and 336 dynamic genes in the brain (2.2%), heart (6.4%), liver (18.3%), muscle (5.5%), and pancreas (1.8%), respectively (Figure 1B). These results highlight the pancreas and the liver with the lowest and the highest number of genes changing their expression with age, respectively. When considering the distribution of the breakpoints of each tissue’s top dynamic genes, we found that, for all tissues except for the muscle, the majority of the changes in gene expression throughout the lifespan occur between mature adulthood and middle age. Until 15 months (approximately the end of middle-age , we observed 64.5% of breakpoints in the brain, 68.4% in the heart, 65.1% in the liver, 41.9% in the muscle, and 79.1% in the pancreas (Figure 1C). Moreover, the median breakpoint times for the brain (12), heart (14), liver (12) and pancreas (15) were not significantly different between these tissues but differ significantly from the median breakpoint time of the muscle (18) (Figure 1C; Supplementary Table 1). Nevertheless, with the exception of the pancreas, the majority of Trendy genes exhibit monotonic expression behaviors (no breakpoints, k=0; Figure 1D), either continuously increasing or decreasing expression throughout time. Biotype assessment of the top dynamic genes in each tissue based on the Ensembl classification (see Methods – Data pre-processing and normalization) showed that dysregulation mainly occurs at the protein coding level (Figure 1E; Supplementary Table 2). The observed under-enrichment of long non-coding RNAs (lncRNAs) in each tissue when compared with the reference genome, might be an artifact resulting from the filtering process as these RNA genes tend to exhibit lower expression levels (Supplementary Table 2). Nonetheless, changes in the expression of long non-coding RNAs may be of relevance, especially in the brain, liver, muscle, and pancreas (Figure 1E).
Different subsets of co-expressed genes exhibit specific age- and sex-related trajectories
As a second approach to defining age-regulated genes, we performed weighted gene correlation network analysis (WGCNA) to explore age- and sex-associated co-expression patterns of gene expression. Similar to the regression analysis, WGCNA was performed on each tissue independently, with each tissue’s co-expression network comprising a variable number of modules of positively correlated genes (Supplementary Table 3).
One of the great advantages of WGCNA is that this methodology takes into consideration the inter-dependency of gene expression. Rather than analyzing each gene independently, WGCNA groups together correlated genes into modules, taking advantage of a power transformation on the pairwise gene correlations to accentuate strong correlations and play down weak ones. The resulting co-expression modules can then be associated with traits of interest (such as age or sex) and the strength of the association can be used to select modules of interest. WGCNA is also very useful for gene prioritization as it allows to rank genes within each cluster and identify intramodular hub genes, likely to be important for the phenotype .
To select the most interesting modules for the aging process, as well as sex-dimorphic genes, we followed a two-step approach. First, for each module, the corresponding gene expression profiles were summarized into a representative - module eigengene (ME; see Methods - Identification of significantly age- and sex- associated modules, hub genes, and Trendy-module-hub overlapping genes) - and this illustrative profile was correlated with age and sex. All modules with a significant (false discovery rate (FDR) < 0.05) and at least moderate (≥ 0.5) correlations between their ME and age were selected. The brain and the heart both exhibited 3 modules significantly associated with age (2 positive and 1 negative in the brain; 1 positive and 2 negative in the heart); the liver displayed 4 age-associated modules (3 positive and 1 negative); the muscle exhibited 5 modules correlated with age (3 positive and 2 negative); and the pancreas did not have any modules significantly correlated with age (Figure 2A). Regarding association with sex, the brain exhibited only 1 positively sex-correlated module; the muscle showed 4 sex-associated modules (2 positive and 2 negative) while the heart and the pancreas exhibited none (Figure 2A). In line with what we previously observed, the liver appears to be the analyzed tissue whose gene expression is the most influenced by sex, as it displays 8 modules significantly sex-associated (4 positive and 4 negative) (Figure 2A; Supplementary Figure 1). Because the pancreas did not exhibit any age- or sex- associated modules, it was not considered for further analyses (Figure 2A). Noteworthily, none of the age-associated modules overlapped with sex-associated ones, with the exception of the Darkgrey module in the liver (Figure 2A). Because sex was treated as a binary variable, with 0 encoding females and 1 encoding males, all positive associations in sex-associated modules correspond to over-representation in males, while all negative-associations correspond to over-representation in females.
Figure 2. Weighted gene co-expression network module selection. (A) Correlation between each module’s eigengene (ME) with age and sex. Each tissue exhibits a variable number of modules of co-expressed genes (brain: 24; heart: 19; liver: 36; muscle: 21; pancreas: 10), and unassigned genes are clustered together in the grey module (not shown). ME is the first principal component of the expression matrix of a module, thus being the most representative gene expression profile of that group of correlated genes. Cells are annotated with bicor values and corresponding FDR adjusted p-values (inside brackets). Red and blue cells depict positive and negative correlations, respectively. The intensity of color represents the degree of correlation. All modules whose ME’s correlation with the trait of interest is significantly equal or higher than 0.5 were considered (moderate correlation and above; FDR < 0.05; significant correlations with age marked with * and significant correlations with sex marked with º). (B) Correlation between module membership (MM) and gene significance (GS) of the previously selected significant modules. MM is obtained by correlating the expression of individual genes to the ME, and GS corresponds to the absolute value of the correlation between individual genes and the trait of interest. Only modules with moderate or higher (≥ 0.5) and significant (p-value < 0.05) correlations with age and sex were considered for subsequent analysis (marked with * and º, respectively).
Next, the gene significance (GS) and module membership (MM) measures were analyzed. GS refers to the absolute value of individual correlations of genes to the trait of interest, whereas MM relates to the individual correlations of genes to the ME. High correlations between these two measures are indicative of genes that are highly significant for the aging process and for sex differences in gene expression, being as well highly important to the module. Modules exhibiting significant (p-value < 0.05) and at least moderate MM-GS correlations (≥ 0.5) were selected. From the 3 previously selected age-associated modules in the brain, only 1 exhibited significant MM-GS correlations Tan module), and the single selected sex-associated module also displayed a significant MM-GS correlation (Grey60 module) (Figure 2B). In the heart, 2 out of 3 age-associated modules displayed significant MM-GS correlations (Tan and Blue modules; Figure 2B). As for the liver, from the 4 previously selected age-related modules, 2 exhibited significant MM-GS correlations Salmon and Darkturquoise modules), while from the 8 selected sex-associated modules, 7 also passed the MM-GS correlation criteria Tan, Red, Darkolivegreen, Darkgrey, Cyan, Brown and Blue modules) (Figure 2B). Interestingly, the genes in the Darkgrey module did not seem to be highly associated with aging as the MM-GS correlation is rather small (0.21; Figure 2B). For this reason, we will consider the Darkgrey module as only a sex-associated module for further analyses. Finally, in the muscle, 2 out 5 age-related modules were selected based on MM-GS correlations (Magenta and Brown modules), and all 4 sex-associated modules also exhibited significant MM-GS correlations (Red, Purple, Greenyellow and Blue modules) (Figure 2B). These selected clusters of genes can be considered the most relevant to the aging process and for sex-dimorphic gene expression.
In order to evaluate each module’s age- and sex-related gene expression, we plotted a heatmap of VST-normalized expression values of the genes belonging to the module, accompanied by a bar plot representing the ME expression profile (Figure 3A; Supplementary Figure 2). This visualization allowed us to better understand the behavior of the age- and sex-related genes across time, as well as to identify the lifespan periods where shifts in expression occur. The brain showed an increasing trend in gene expression with age (Tan module; bicor=0.86; p-value=3e-13), with the shift from down- to up-regulation occurring around the transition from middle- to old-age (15 to 18 months) (Figure 3A, upper left panel). Regarding the genes present in its sex-associated module – Grey60 (bicor=0.84; p-value=9e-12), they are over-expressed in males, and no obvious trend of increased/decreased expression throughout time was observed (Supplementary Figure 2, upper left panel). As for the heart, the 2 significant age-related modules display different expression trends. The Blue module (bicor=-0.57; p-value=5e-04) exhibits a decreasing trend in gene expression with age, and the transition from up- to down-regulation happened at old-age (18 to 21 months) (Figure 3A, upper left panel). Conversely, the Tan module presents an increasing trend in expression throughout the lifespan, transitioning from down- to up-regulation within middle-age (12 to 15 months) (Figure 3A, upper central panel). Concerning the liver, both age-associated modules displayed an increasing trend in gene expression over time, the difference being the onset of expression change. In the Salmon module (bicor=0.71; p-value=1e-06), gene expression starts to increase within middle-age (9 to 12 months), whereas in the Darkturquoise (bicor=0.54; p-value=0.001) this shift is less defined, probably occurring around the transition from mature adulthood to middle-age (6 to 9 months) (Figure 3A, lower left panels). In the case of modules significantly correlated with sex in the liver, the Brown (bicor=-0.89; p-value=1e-14), the Darkolivegreen (bicor=-0.53; p-value=7e-04) and the Red (bicor=-0.68; p-value=2e-06) modules are female-associated, and none displayed any obvious trends in gene expression changes over time (Supplementary Figure 2, upper right and central panels). As for the male-associated modules in the liver, all exhibit slightly decreasing trends in the expression of their genes, with all shift points in expression occurring within old age. Both the Blue (bicor=0.84; p-value=1e-11) module and the Darkgrey (bicor=0.79; p-value=1e-09) modules exhibited a shift between 24 and 27 months (within old age) (Supplementary Figure 2, upper right and central panels). Moreover, the Cyan module (bicor=0.55; p-value=6e-04) displayed its shift in gene expression around the transition from 18 to 21 months (Supplementary Figure 2, upper right and central panels), whereas the Tan module’s (bicor=0.55; p-value=6e-04) expression change occurred between 21 and 24 months (Supplementary Figure 2, upper right and central panels). Concerning the muscle, we observed opposite trends in the two age-associated modules. In the Magenta module (bicor=0.73; p-value=3e-07) we observed increased gene expression across the lifespan, whereas in the Brown module (bicor=-0.76; p-value=8e-08), the expression tends to decrease over time (Figure 3A, lower right panels). Interestingly, in both modules the shift in expression occurs within middle-age, around the transition from 9 to 12 months (Figure 3A, lower right panels). Finally, in regards to the sex-associated modules in the muscle, only the Blue module (bicor=0.5; p-value=0.004), which is male-related, displayed a discernible trend of increased expression across the lifespan, with the transition from down- to up-regulation taking place within middle-age (9 to 12 months) (Supplementary Figure 2, lower left panel). The other 3 sex-related modules in this tissue did not exhibit any obvious trends in gene expression over time, with the Red module (bicor=0.86; p-value=3e-12) being male-associated, and the Greenyellow (bicor=0.85; p-value=4e-12) and Purple (bicor=0.76; p-value=2e-08) modules being female-associated (Supplementary Figure 2, lower right panels).
Figure 3. Weighted gene co-expression network significantly age-associated modules. (A) Gene expression profile of each significantly age-associated module. The heatmaps (top) display the standardized expression (z-score) of individual genes (rows) per sample (columns), whereas the bar plots (below) represent the ME expression profile. Each bar of the bar plot corresponds to the same samples of the heatmap. Negative (positive) values of ME expression relate to the under-expression (over-expression) of genes in each module’s heatmap (blue and red colors, respectively). (B) Intramodular hub gene identification. For each module, genes with individual GS > 0.2 and MM > 0.8 were considered to be the most functionally important (inside grey rectangles).
Furthermore, in each module, we identified the genes with the highest MM and GS – Hub genes - as they are important elements of the module, as well as the most significantly associated with the trait of interest (see Methods - Identification of significantly age- and sex-associated modules, hub genes, and Trendy-module-hub overlapping genes; Figure 3B and Supplementary Figure 3; and Supplementary Table 4).
Altered genes and biological networks provide tissue-specific markers of aging
In order to integrate the results from the two described approaches for finding age-dysregulated genes, we intersected the resulting gene lists and evaluated their functional implications. For each tissue, we compared the previously identified Trendy top dynamic genes with the gene sets of interest resulting from WGCNA (Module and Hub genes) and selected for functional analysis the intersection of the Trendy and the Hub genes (Figure 4; Supplementary Figure 4; Supplementary Table 5). As expected, all Hub genes overlapped with the Module genes. Because our interest is to explore gene expression dysregulation across the lifespan, we performed this analysis on all age-associated modules and on only the sex-associated modules exhibiting discernible trends of increased/decreased gene expression over time (i.e., the Blue, Cyan, Darkgrey, and Tan modules in the liver, and the Blue module in the muscle).
Figure 4. Gene overlap between Trendy genes, module genes, and hub genes in the age-associated modules with significantly enriched GO terms. Bars represent intersection size and colored circles depict the gene sets involved. Genes in common in the Trendy and hub gene sets were considered for further analysis (identified with *).
In the brain Tan module, 34 genes were selected for further analysis (Figure 4; Supplementary Table 5). As for the heart, the size of the overlap between Trendy and Hub genes was 11 and 94, for the Tan and Blue modules, respectively (Figure 4; Supplementary Table 5). In the liver, regarding age-associated, 5 and 9 intersecting genes were selected from the Salmon and Darkturquoise modules, respectively (Figure 4; Supplementary Table 5), whereas 92, 12, 20 and 13 genes were selected from the Blue, Cyan, Darkgrey, and Tan sex-associated modules, respectively (Supplementary Figure 4; Supplementary Table 5). Finally, in the muscle we proceeded the analysis with 9, 31, and 21 genes from the Magenta, Brown and Blue modules, respectively (Figure 4 and Supplementary Figure 4; Supplementary Table 5).
To better understand the functions underlying these signature gene lists, we performed an enrichment analysis on GO BPs and selected the ones with an FDR adjusted p-value of less than 0.05 (see Methods - Functional characterization of Trendy-module-hub genes’ overlap; Supplementary Table 6). To deal with the large number of significant BPs exhibited by some modules and to provide a clear picture of how they and their associated genes relate to each other, we constructed a network of these results, with nodes representing GO terms and edges corresponding to the overlap of genes between any two nodes (see Methods - Network visualization of functionally enriched terms). After having constructed the network, we addressed GO term redundancy by clustering together nodes based on gene overlap similarity and then assigning automatically created labels from the most frequent words in the cluster, as well as words adjacent to the most represented ones (see Methods - Network visualization of functionally enriched terms; Supplementary Table 6) [as seen in 53, and in 54]. To further simplify the visualization, we created a summary network based on the generated clusters, where all the nodes belonging to the same cluster collapsed into a meta-node and all the edges connecting the different clusters collapsed into meta-edges as well  (Supplementary Table 6). Selected results are present in Table 1.
Table 1. Tissue-specific age-dysregulated functions.
|Tissue||Module||Summary network meta-node||Summary network FDR||No. of nodes||Meta-node choice criteria||Genes involved||WGCNA expression trend||WGCNA expression shift|
|Brain||Tan||mediated immunoglobulin regulation response||8.06E-03||142||Largest||B2m, C1qa, C1qb, C1qc, C3, C4b, Csf1, Ctsh, Ctss, Ctsz, Cx3cr1, Gfap, H2-D1, H2-K1, H2-T23, Hexb, Il33, Irf7, Itgb2, Lag3, Lgals3, Neat1, Psmb8, Serpina3n, Slc11a1, Tap2, Tapbp||Increasing with age||Middle to old-age|
|Brain||Tan||neuron death activation microglial||6.59E-03||7||Most significant||C1qa, Csf1, Ctsz, Cx3cr1, Il33, Slc11a1||Increasing with age||Middle to old-age|
|Heart||Blue||branched-chain amino process metabolic||1.57E-02||17||Largest||Acaa2, Adk, Auh, Ckm, Ckmt2, Cyc1, Eno3, Gnpat, Hmgcl, Isca1, Mccc2, Nfs1, Pdhb, Suclg1||Decreasing with age||Old-age|
|Heart||Blue||acid catabolic lipid fatty||4.19E-04||12||Most significant||Acaa1a, Acaa2, Adipor1, Auh, Eci2, Etfb, Gnpat, Hadh, Hmgcl, Mccc2, Pgm2, Phyh, Ptges2, Smpd1||Decreasing with age||Old-age|
|Liver||Darkturquoise||response immune regulation cell||7.07E-03||57||Both||Cd19, Cd79a, Cd79b, Ighm, Igkc, Igkv3-5, Iglc2, Jchain, Mzb1||Increasing with age||Adulthood to middle-age|
|Liver||Salmon||regulation assembly cell positive||3.32E-02||184||Both||Ccl5, H2-Aa, H2-Eb1, Ntrk2, Slamf7||Increasing with age||Middle-age|
|Liver||Blue||lipid biosynthetic process||3.43E-02||1||Both||Cdipt, Elovl2, Elovl3, Hsd17b12, Mecr, Pip5k1a, Pisd, Rest, Scp2, Serinc1, Smpd1||Decreasing with age (males)||Old-age|
|Liver||Cyan||endoplasmic reticulum golgi vesicle-mediated||2.80E-03||4||Most significant||Arcn1, Copg1, Sec22b, Sec24d||Decreasing with age (males)||Old-age|
|Liver||Cyan||protein transmembrane response development||1.41E-02||21||Largest||Arcn1, Copg1, Hspa5, Sdf2l1, Sec22b, Sec24d, Sec61a1, Serp1||Decreasing with age (males)||Old-age|
|Muscle||Brown||regulation morphogenesis negative development||1.12E-02||82||Both||Angptl1, Anxa2, Cd34, Col3a1, Col5a1, Col5a2, Col6a1, Col6a2, Col6a3, Dok2, Fn1, Igfbp6, Itgbl1, Lrp1, Ndn, Nid1, Pi16, Serpinf1, Serping1, Sparc, Ssc5d, Tgfbi, Timp2||Decreasing with age||Middle-age|
|Muscle||Blue||proteasomal ubiquitin-independent protein catabolic||1.78E-02||1||Most significant||Psmb1, Psmb3||Increasing with age (males)||Middle-age|
|Muscle||Blue||ribonucleoside triphosphate metabolic process||2.89E-02||13||Largest||Atp5l, Atp5g3, Cox5a, Ndufa8||Increasing with age (males)||Middle-age|
|The choice of the representative meta-nodes was based on the number of nodes (GO terms) included (Largest) and/or on the p-values (Most significant). In the sex-associated modules, the sex in which the genes are enriched is indicated inside parentheses. Relates to Supplementary Table 6.|
Neuronal death and immune response processes are upregulated in the aging brain
GO enrichment analysis of the overlap of the brain Trendy genes with the Tan module hub genes identified 233 significantly over-represented BPs, allocated into 9 meta-nodes (Supplementary Table 6). The largest meta-node is involved in antigen-mediated immunity and comprises 142 GO terms, thus representing more than half of all significant processes observed for this set (Brain Tan module: Table 1; Supplementary Table 6). The genes present in this group exhibit increased expression with increasing age, shifting from down- to up-regulation around the transition from middle to old age (15-18 months; Brain Tan module: Figure 3A; Table 1).
Furthermore, 7 processes related to neuronal death were found to be the most significantly enriched in the aging brain (Brain Tan module: Table 1; Supplementary Table 6). This meta-node also includes genes whose expression tend to increase throughout aging with the shift from down- to up-regulation occurring from middle- to old-age (15-18 months; Brain Tan module: Figure 3A; Table 1).
Decline of cardiac ribonucleotide and fatty acid metabolism with age
In the heart, we identified 11 clusters of similar GO terms, comprising 79 significantly enriched GO BPs based on the intersection between the heart top dynamic genes and the heart Blue module hub genes (Heart Blue module: Table 1; Supplementary Table 6). The most significant meta->node includes processes related to fatty acid metabolism, whereas the largest cluster comprises BPs related to ribonucleotide metabolism (Heart Blue module: Table 1; Supplementary Table 6). The genes found to be involved in these processes exhibit decreased expression along the lifespan, with the transition from up- to down-regulation occurring late in life (18-21 months; Heart Blue module: Figure 3A; Table 1).
Muscle regeneration is impaired with aging and mitochondrial energy metabolism and proteasomal activity are altered in the muscle of aged male mice
Over-representation analysis of GO terms in the gene set obtained from the overlap of the muscle Trendy genes and the muscle Brown module hub genes resulted in a highly interconnected network of 82 BPs organized into a single meta-node (Supplementary Table 6). This cluster includes genes related to muscle tissue regeneration (Muscle Brown module: Table 1; Supplementary Table 6), which are down-regulated across the lifespan. A shift from up- to down-regulation occurs in middle age (9-12 months; Muscle Brown module: Figure 3A; Table 1).
Additionally, the muscle exhibits one gene module whose expression is male-associated where we could observe an increasing trend in expression over time (Muscle Blue module: Supplementary Figure 2; Table 1). After performing GO over-representation analysis on the overlap of muscle dynamic genes with the muscle Blue module hub genes we identified 28 processes allocated into 4 meta-nodes (Muscle Blue module: Supplementary Table 6). The largest cluster is related to the metabolism of adenosine triphosphate (ATP), while the most significant one is associated with proteasomal activity (Muscle Blue module: Table 1; Supplementary Table 6). The genes involved in these processes exhibit increased expression in aged males, with the transition from under- to over-expression taking place around middle-age (9-12 months; Muscle Blue module: Supplementary Figure 2; Table 1).
Liver aging is characterized by a global dysregulation of immune function and sex-specific differences in lipid metabolism and stress response activity
In the liver, a total of 247 BPs were identified as being significantly enriched among the gene lists resulting from the overlap of the liver top dynamic genes and the age-associated modules hub genes (Liver Darkturquoise and Salmon modules: Table 1; Supplementary Table 6). In the Trendy- Darkturquoise module gene set, 63 enriched GO terms were allocated into 2 meta-nodes, with the largest and most significant (n = 57) relating to immune cell activation, differentiation and proliferation (Liver Darkturquoise module: Table 1; Supplementary Table 6). The genes involved in these processes exhibit increased expression throughout the lifespan, with the shift from down- to up-regulation occurring between adulthood to middle age (6-9 months; Liver Darkturquoise module: Figure 3A; Table 1). Regarding the Trendy-Salmon module hub gene overlap, 1 cluster of processes also related to immune response was identified, comprising 184 GO terms (Liver Salmon module: Table 1; Supplementary Table 6). The identified genes within this cluster present an increased trend in expression with increasing age, shifting from down- to up-regulation within middle age (9-12 months; Liver Salmon module: Figure 3A; Table 1).
Similar to what we have observed in the muscle, the liver also exhibits sex-dimorphism in gene expression with aging. In the liver Blue module, GO over-representation analysis on the overlap of top dynamic genes with hub genes allowed for the identification of 1 BP related to biosynthesis of lipids, whose comprising genes exhibit decreased expression over time, shifting from up- to down-regulation really late in life, around the transition from 24 to 27-months (Liver Blue module: Supplementary Figure 2; Table 1; Supplementary Table 6). Moreover, in the Cyan module, 25 identified processes were allocated into 2 meta-nodes, the most significant one related to the endoplasmic reticulum (ER) transport system, and the largest one associated with ER stress response (Liver Cyan module: Table 1; Supplementary Table 6). The genes involved in these processes exhibit decreased expression in aged males, with the transition from over- to under-expression taking place within old-age (18-21 months; Liver Cyan module: Supplementary Figure 2; Table 1).
Similar biological processes altered across multiple tissues
Evidence from other studies have showed that differences in gene expression across the lifespan are strongly associated with tissue type. For that reason, all the analyses were performed in each tissue independently. Our data confirmed this trend since few players that are associated with aging are shared between tissues (5 genes in common between only the heart and the liver) (Figure 5 - left panel; Supplementary Table 7). Nonetheless, when we compared the GO terms enriched between tissues, we found a much higher overlap than that observed at the gene level (Figure 5 - right panel; Supplementary Table 8). The highest number of shared BPs (43) was observed between the brain and the liver, followed by 6 processes shared between the heart and the muscle. The brain and the muscle also exhibited 6 GO terms in common, while the heart and the liver shared 3 BPs. Finally, the muscle and the liver, presented only 1 BP in common (Figure 5 – right panel; Supplementary Table 8). As described above, we organized the GO terms in summary networks (Supplementary Table 8) and summarized the results in Table 2.
Figure 5. Overlap of genes and biological processes between the brain, heart, liver, and muscle. Upset plots depicting the gene overlap between Trendy genes, module genes, and hub genes per tissue (left), as well as the overlap of the enriched GO terms in the same tissues (right). Bars represent intersection size and colored circles depict the gene/GO term sets involved. Each tissue's gene list results from the intersection of Trendy, module, and hub genes. For this plot we considered all significantly age-associated modules and the significantly-sex associated modules with an observable increase/decrease in gene expression over time (i.e., liver tan, blue, cyan and Darkgrey modules, and muscle blue module). In tissues with more than one significant module (i.e. the heart, the muscle and the liver), the gene list results from the combination of each module's intersection, and the GO term list results from the combination of each module's GO terms. GO terms in common at least in two tissues were considered for further analysis (identified with *).
Table 2. Inter-tissue age-dysregulated functions.
|Related function||Summary network description||No. of nodes||Genes involved - Brain||Genes involved - Liver||WGCNA expression trend - Brain||WGCNA expression trend - Liver||WGCNA expression shift - Brain||WGCNA expression shift - Liver|
|Immune response||Response immune mediated immunoglobulin||13||B2m, C1qa, C1qb, C1qc, C3, C4b, Ctsh, Gbp3, H2-D1, H2-K1, H2-T23, Il33, Irf7, Itgb2, Lag3, Lgals3, Lyz2, Slc11a1, Tap2||Cd19, Cd79a, Cd79b, H2-Aa, H2-Eb1, Ighm, Igkc, Iglc2, Jchain, Slamf7||Increase||Increase||Middle to old-age|
|Adulthood to Middle-age|
6-9 months and
|regulation production cell proliferation||26||B2m, C1qa, C1qc, C3, Csf1, Ctsz, Cx3cr1, H2-D1, H2-K1, H2-T23, Il33, Itgb2, Lag3, Lgals3, Slc11a1, Tap2||Ccl5, Cd19, Cd79a, H2-Aa, Ighm, Igkv3-5, Mzb1, Ntrk2, Slamf7||Increase||Increase||Middle to old-age|
|processing presentation peptide antigen||4||B2m, Ctss, H2-D1, H2-K1, H2-T23, Psmb8, Slc11a1, Tap2, Tapbp||H2-Eb1, H2-Aa||Increase||Increase||Middle to old-age|
|Related function||Summary network description||No. of nodes||Genes involved - Brain||Genes involved - Muscle||WGCNA expression trend - Brain||WGCNA expression trend - Muscle||WGCNA expression shift - Brain||WGCNA expression shift - Muscle|
|Tissue regeneration||angiogenesis||1||C3, Lgals3, Itgb2, Ctsh, Cx3cr1||Fn1, Serpinf1, Sparc, Cd34, Anxa2||Increase||Decrease||Middle to old-age|
|glial cell migration||1||Csf1, Hexb||Fn1, Ndn, Lrp1||Increase||Decrease||Middle to old-age|
|extracellular structure organization response||3||C3, Ctss, Gfap, Lgals3, Neat1, Slc11a1||Anxa2, Cd34, Col3a1, Col5a1, Col5a2, Fn1, Lrp1, Nid1, Serping1, Tgfbi||Increase||Decrease||Middle to old-age|
|regulation peptidase activity||1||Serpina3n, Ctsh, Psmb8, Ctsd||Serpinf1, Pi16, Serping1, Timp2||Increase||Decrease||Middle to old-age|
|Related function||Summary network description||No. of nodes||Genes involved - Heart||Genes involved - Muscle||WGCNA expression trend – Heart||WGCNA expression trend - Muscle||WGCNA expression shift - Heart||WGCNA expression shift - Muscle|
|Energy metabolism||respiratory electron transport atp||6||Coq9, Cyc1, Eno3, Ndufa10, Ndufs2, Sdhd, Slc25a12, Uqcrc1||Atp5l, Atp5g3, Cox5a, Ndufa8||Decrease||Increase (males)||Old-age|
|Related function||Summary network description||No. of nodes||Genes involved - Liver||Genes involved - Heart||WGCNA expression trend – Liver||WGCNA expression trend - Heart||WGCNA expression shift - Liver||WGCNA expression shift - Heart|
|Glucocorticoid signaling||receptor signaling pathway glucocorticoid||3||Ntrk2||Ppp5c, Phb||Increase||Decrease||Middle-age|
|Related function||Summary network description||No. of nodes||Genes involved - Liver||Genes involved - Muscle||WGCNA expression trend – Liver||WGCNA expression trend - Muscle||WGCNA expression shift - Liver||WGCNA expression shift - Muscle|
|Response to amino acid stimulus||cellular response to amino acid stimulus||1||Ntrk2||Col3a1, Col5a2, Col6a1||Increase||Decrease||Middle-age|
|In the sex-associated modules, the sex in which the genes are enriched is indicated inside parentheses. Relates to Supplementary Table 8.|
The liver and the brain share 43 GO biological processes clustered into 3 meta-nodes (Brain-Liver; Table 2; Supplementary Table 8). In line with what we observed at the tissue-level, immune processes are represented in the intersection of GO terms between these two tissues. In both tissues, the genes involved in the processes comprising these meta-nodes exhibit an increasing trend of expression across the lifespan (Brain Tan and Liver Salmon and Darkturquoise modules: Figure 3A; Brain-Liver: Table 2; Supplementary Table 8). However, the dysregulation of the genes involved in these processes occurs later in the brain, in the transition from middle to old age (15-18 months), than in the liver (6-9 months and 9-12 months; Brain Tan and Liver Salmon and Darkturquoise modules: Figure 3A; Brain-Liver: Table 2; Supplementary Table 8).
Interestingly, in the brain and the muscle the commonly affected BPs are related to tissue regeneration, with the 6 identified GO terms allocated to 4 clusters (Brain-Muscle: Table 2; Supplementary Table 8). In the brain, the genes comprised by these meta-nodes display increased expression across the lifespan, whereas the corresponding muscle genes displayed a down-regulation trend (Brain Tan and Muscle Brown modules: Figure 3A; Brain-Muscle: Table 2). Moreover, the shifts in regulation occurred late in life in the brain, around the transition from middle to old age, whereas in the muscle they are observed earlier, within middle age (15-18 and 9-12 months, respectively; Brain Tan and Muscle Brown modules: Figure 3A; Brain-Muscle: Table 2).
In the case of the heart and muscle, the 6 age-related dysregulated processes are related to energy metabolism, particularly to the synthesis of ATP (Heart-Muscle: Table 2; Supplementary Table 8). Notably, the overlap of these processes is observed between the hearts of both male and female mice, and the muscle of only male mice. In the heart, the genes related to these processes exhibit a decreasing trend in expression, transitioning from over- to under-expression within old age (18-21 months; Heart Blue module: Figure 3A; Heart-Muscle: Table 2;Supplementary Table 8). As for the genes in the muscle they tend to increase their expression in males over time, shifting from down- to up-regulation within middle age (9-12 months; Muscle Blue module: Supplementary Figure 2; Heart-Muscle: Table 2; Supplementary Table 8).
The liver and the heart share 3 biological processes, all related to glucocorticoid signaling (Liver-Heart: Table 2; Supplementary Table 8). Interestingly, the genes enriched in these GO terms display opposite trends in expression throughout life. In the liver, we have observed increased expression with age, with a shift point within middle age (9-12 months; – Liver Salmon module: Figure 3A; Liver-Heart: Table 2; Supplementary Table 8), whereas in the heart, we found gene expression to decrease across the lifespan, transitioning from up- to down-regulation later in life (18-21 months; Heart Blue module: Figure 3A; Liver-Heart: Table 2; Supplementary Table 8).
As for the liver and the muscle, we found only 1 process in common between these tissues involving the response to amino acid stimulus (Liver-Muscle: Table 2; Supplementary Table 8). Despite the genes comprised in this single meta-node changing their expression within middle age (9-12 months) in both tissues, they present opposite expression trends, increasing in the liver and decreasing in the muscle (Liver Salmon and Muscle Brown modules: Figure 3A; Liver-Muscle: Table 2; Supplementary Table 8).
In this work, we provide an in-depth characterization of the age-associated alterations in gene expression throughout the murine lifespan, through the re-analysis of publicly available mouse aging transcriptomic data (GSE132040). We have combined, for each tissue, a segmented regression model fitted to each gene to select the top dynamic genes over time with a network-based approach to establish tissue-specific clusters of co-expressed genes correlated with increasing age and to identify sex-dimorphic gene expression patterns. By integrating the results of these two methodologies, we were able to both prioritize aging genes and establish gene expression signatures of aging with greater confidence. Among the selected tissues, the liver showed the highest percentage of dynamic genes across the lifespan, which is in line with recent works from other laboratories  and may be related to its high metabolic activity and protein synthesis rate [56, 57]. Nevertheless, these observations do not consider the magnitude of the dysregulation so conclusions must be drawn carefully as only one aspect of age-related gene expression dysregulation was addressed. We have also observed that most alterations in gene expression occur between 12 and 15 months in all tissues but the muscle, where the majority of changes took place in old age. Interestingly, Schaum and colleagues also observed a late-life gene expression dysregulation in this tissue, with the largest number of differentially expressed genes (DEGs) occurring in the 24 and 27 month age groups . Additionally, the pancreas appeared to be largely unaffected by aging since it displayed the lowest amplitude of dysregulation and no significantly age-associated gene co-expression modules, also corroborating the findings of Schaum et al. . The reason for this is not clear and data analysis does not highlight technical issues with the pancreas RNA-Seq dataset, however, these results should be interpreted with caution since RNA isolation may be affected by the high level of ribonucleases present in this tissue [58, 59].
The observed gene expression alterations across the lifespan generally reflect loss of tissue function and homeostasis (Supplementary Discussion), with different tissues exhibiting diverse onsets of gene expression dysregulation, as has been described before [6, 20]. We have also observed the existence of significant sex-dimorphic expression patterns in the brain, muscle, and liver, with the former exhibiting the lowest number of sex-biased modules, and the latter displaying the highest. This observation is concordant with previous research in mice that showed the highest dimorphic gene expression to occur in the liver, followed by the muscle, and then by the brain . These findings are also in line with those from the original study of the dataset re-analyzed in this work, which described prominent sex effects in gene expression in the liver and other tissues, not mentioning, however, neither the brain nor the muscle . Sex-differences in gene expression have already been reported for the heart  and the pancreas  in healthy humans and mice, respectively, however, we did not find modules significantly correlated with sex in the heart and in the pancreas. One limitation of our study, which is inherent to the dataset, is the unbalanced number of subjects from both sexes that may mask more subtle, but existing, sex-differences in gene expression in these tissues.
Genes involved in proteostasis-related pathways were significantly enriched in males over females and exhibited expression differences over time in the liver and the muscle of only male mice. For instance, in the liver, we observed sex-dimorphism in the expression of genes related to ER stress response, with decreased expression in males within old age (18-21 months). Notably, a recent study focusing on sex differences in hepatic ER stress described an increase in markers of the unfolded protein response (UPR) in the liver of male rats in the transition from the prepubertal phase to adulthood . In line with these findings, we observed a peak in the expression of ER-stress related genes in the adult age prior to the described decrease throughout the lifespan. For that reason, further studies addressing sex-differences in these markers in a larger number of time points would be very important for deepening the knowledge on age-related ER-stress. As for the muscle, genes encoding subunits of the proteasome were found to be enriched in males and exhibited increased expression with aging, shifting from under- to over-expression around middle-age (9-12 months). A recent study has reported females as having significantly higher proteasomal activity in several tissues , however, this study did not include muscle tissue.
One major implication from our results is that sex-specific alteration in genes involved in proteostasis-related processes can be detrimental to male survival and longevity since these pathways maintain protein quality control and prevent aberrant protein aggregation and proteotoxic stress [recently reviewed by us in 64]. It is well known that proteostasis decline is a hallmark of aging , leading to proteome imbalances and contributing to protein aggregation, including of amyloid-like aggregates [66–69], particularly in the muscle [70, 71]. Nonetheless, despite past studies having reported age-related declines in proteostasis in mammals, sex-dimorphism in these processes across the lifespan has been less explored, making this a very interesting subject to address in future research.
Additionally, during liver aging there was a decreased male-enriched expression of genes involved in the biosynthesis of lipids (shift point within old age; 24-27 months) while the muscle showed increased male-biased expression (shift point within middle age; 9-12 months) of genes related to energy metabolism, particularly to the synthesis of ATP (Supplementary Discussion).
Although very few commonalities in age-related dysregulation were observed at the gene-level, we realized that the studied tissues could be grouped according to shared biological processes affected by aging. For the most part, the observed shifts in gene expression occurred within middle-age (9-12 months), as seen in the muscle and liver, however, in the brain and the heart these shifts occurred later in life (15-18 months and 18-21 months, respectively). The brain and the liver had common age-related activation of immune responses, with increased expression of immune-related genes occurring later in the brain (15-18 months) than in the liver (6-9 months and 9-12 months). This activation pattern is consistent with the findings of Schaum et al. , while also being implicated in age-related neurodegeneration, as well as in several hepatic age-related diseases, such as non-alcoholic fatty liver disease (NAFLD) and hepatocellular carcinoma (HCC) (Supplementary Discussion). Interestingly, a recent study comparing age-related transcriptomic alterations between different tissue-resident macrophages has shown that the most prominent dysregulation occurred in microglia, followed by Kupffer Cells (KC), the brain and the liver’s resident macrophages, respectively [72; preprint].
Another major implication from our observations is that, similarly to the shift points, the direction of dysregulation varies between tissues, despite shared biological pathways. For instance, the brain exhibits an up-regulation of genes involved in regeneration processes at old age, while the muscle shows an earlier (middle-age) down-regulation of these processes. Interestingly, serpin family genes are altered in both tissues, with the downregulation of Serpinf1 (serine (or cysteine) peptidase inhibitor, clade F, member 1) and Serping1 (serine (or cysteine) peptidase inhibitor, clade F, member 1) in the muscle contrasting with the up-regulation of Serpina3n (serine (or cysteine) peptidase inhibitor, clade A, member 3N) in the brain. Serpins are a family of serine (or cysteine) protease inhibitors involved in several biological functions including homeostasis control . In the muscle, both Serpinf1 and Serping1 are involved in muscle growth and function through regulation of Akt and FoxO signaling pathways [74–76] whereas Serpina3n, which has been linked with increased immune response activity, is upregulated in aging astrocytes throughout the brain , as well as in a Prion disease mouse model , and has been implicated in Alzheimer’s disease (AD) [addressed in 79] and in Multiple Sclerosis (MS) disease progression .
Alterations in signaling and cellular response processes were also shared between tissues, with opposing regulation over time in each tissue-pair (Supplementary Discussion). For example, the liver and the heart both showed alterations in genes involved in glucocorticoid signaling, with an early (9-12 months) increase in expression in the liver contrasting with a late life (18-21 months) decrease in expression in the heart. Moreover, the liver and the muscle shared dysregulation of cellular responses to amino acid stimuli with the increase in gene expression in the liver opposing to the decrease in the muscle, albeit the shift point occurring around the same time (9-12 months). Interestingly, the aging heart exhibits a late life (18-21 months) decreased expression of respiratory metabolism genes, while the muscle showed an earlier increased expression of genes involved in the same processes (9-12 months), but only for male mice.
In conclusion, this work opens new research avenues as it highlights many unexplored genes and mechanisms in the context of healthy aging and temporally contextualizes gene expression alterations. We identified tissue-specific key players of aging and addressed the functional implications of their age-related alterations and sex differences in gene expression throughout the lifespan. We also identified groups of tissues based on shared age-affected processes, potentially uncovering tissue axes of common age-related functional dysregulation. We found that proteostasis impairment is a common feature of aging in the muscle and liver of male mice and may dictate sex-differences in lifespan. Because these alterations occur at a transcriptional level and protein abundances can be post-transcriptionally regulated, we are currently working on integrating these results with translatomic and proteomic data to further explore these findings. Additionally, expanding this approach to the other tissues present in the dataset would be pivotal to comprehensively understand this phenomenon and successfully promote healthy aging strategies.
Materials and Methods
The mouse bulk RNA-Seq data used in this study was made publicly available by the Tabula Muris Consortium [6, 45], and is deposited in NCBI’s Gene Expression Omnibus (GEO) under the GEO Series accession number GSE132040 . The original dataset consists of transcriptomic data from 17 male and female mouse tissues across 10 time points (1, 3, 6, 9, 12, 15, 18, 21, 24, and 27 months). For this study, we excluded the 1-month-old samples to avoid the influence of developmental genes , and selected the brain, heart, muscle, liver, and pancreas for further analysis (Supplementary Table 9). The RNA extraction, cDNA library preparation, RNA sequencing, read quality control, pre-processing and alignment, transcriptome reconstruction, and expression quantification steps are reported in the original study, in the GEO webpage (GSE132040 entry), and in the protocols.io repository (2uvgew6 entry) [6, 44, 81]. Briefly, libraries were sequenced on the NovaSeq 6000 Sequencing Systems (Illumina) and originated 100-bp paired-end reads which were then de-multiplexed with bcl2fastq (v.2.20.). Read alignment to the GRCm39.p6 (Mus musculus) reference genome with Gencode (v.M19) annotations was performed using STAR (v. 2.5.2b). Transcript reconstruction and gene expression quantification was performed with HTSeq (v. 0.6.1p1).
Data pre-processing and normalization
Similar to what was described in the original study , we discarded samples with library size smaller than 4.000.000 reads across all genes (Supplementary Table 7). We further carried out outlier sample identification and removal, which was not performed in the original analysis. Outlier samples were identified based on the sample network approach [82, 83] and excluded if their standardized connectivities (z.K) were more than 2 standard deviations away from the mean z.K (Supplementary Table 9; Supplementary Figure 5).
Gene symbols were associated with Ensembl (release 99) biotype annotations  (Mus musculus reference genome GRCm39) using the R package biomaRt (v. 2.44.0) [85, 86]. Low count genes were pre-filtered and only genes with total read count higher than 10 in at least n samples were kept, with n corresponding to each tissue’s minimum group size, whereas in the original study no information regarding gene filtering is provided. Read count data was normalized across samples and transformed with DESeq2’s (v. 1.28.1)  estimateSizeFactors and vst functions , respectively, as described by Schaum and colleagues . PCA based on the 500 genes with highest row-wise variance (i.e., across all samples) was performed to identify the highest contributing sources of variance, but only for the 5 selected tissues. In line with the observations of Schaum et al. , since all samples segregated mainly by tissue (Figure 1A), the subsequent analyses were carried out separately for the brain, heart, liver, muscle, and pancreas (see Methods - Dataset characterization).
Trendy segmented regression analysis
Segmented regression analysis was carried out in normalized expression data (see Methods – Data pre-processing and normalization) using Trendy (v. 1.8.2) . For each gene, a set of segmented regression models were fit, with gene expression as a function of time. We also included ‘Sex’ in the model. Each model comprised a varying number of breakpoints, representing a dynamic change in the gene expression profile, and ranging from 0 (linear regression) to 8 (maxK parameter). We also determined a minimum number of samples present in a segment (expression profile interrupted by a breakpoint), based on the minimum group size in each tissue (brain: 3, heart: 4, liver: 3, muscle: 2, pancreas: 2). Then, the method selected the fitted model having the optimal number of breakpoints by identifying the one with the smallest Bayesian information criterion (BIC) value. We defined top dynamic genes based on the adjusted R2 values (indicative of the goodness of fit of a model) and on the segment slope p-values. The cutoff for the adjusted R2 was defined for each tissue independently (brain: 0.2, heart: 0.2, liver: 0.1, muscle: 0.3, pancreas: 0.3) and chosen as the value above which less than 1% of the genes were kept after a permutation procedure [as performed in 48]. Thus, a gene was considered to be top dynamic (Trendy) if the adjusted R2 was higher than the tissue-specific established cutoff and if at least one segment was significant (segment slope p-value < 0.1). For the comparison of breakpoint distribution between the tissues, we performed a Kruskal-Wallis test, followed by a Dunn’s test of multiple comparisons using the kruskal.test and dunnTest functions, from the R packages stats (v. 4.1.0) and FSA (v. 0.8.32), respectively. For testing biotype enrichment, we performed the Fisher’s Exact Test using the R function chisq.test, from the stats (v. 4.1.0) package.
Gene co-expression network construction and module construction
The WGCNA R package (v. 1.69)  was used to construct co-expression networks for the VST-transformed data (see Methods – Differential gene expression analysis). Due to the large size of the datasets, an automatic block-wise network construction and module detection approach was chosen .
First, genes with zero variance across all samples were flagged and excluded. Then, for each filtered dataset a correlation matrix was calculated based on biweight midcorrelation (bicor) values and raised to specific soft thresholding powers (β; brain: 7; heart: 7; liver: 5; muscle: 6; pancreas: 7; Supplementary Figure 6). In the cases where the scale-free topology fit index failed to reach values above 0.8, the soft-threshold power was chosen based on the number of samples and mean connectivity values . The resultant signed adjacency matrices were used to compute measures of topological overlap between each pair of genes, present in Topological Overlap Matrices (TOM). Next, genes in each dataset were hierarchical clustered (average linkage method) based on topological overlap dissimilarity (1- TOM). Modules of co-expressed genes were constructed accounting for a minimum size of 50 genes, a dendrogram branch merge cut height of 0.15, and default module detection sensitivity (deepsplit = 2) for all datasets except for the liver (deepsplit = 4).
Identification of age- and sex-associated modules, hub genes, and Trendy-module-hub overlapping genes
An initial selection of modules was based on the association of each module eigengene (ME) with age and sex. Age was treated as a continuous variable, while sex was transformed into a binary variable, with 0 encoding females and 1 encoding males. ME is the first principal component of the expression matrix of a module and is usually considered to be the most representative gene expression profile of that group of correlated genes. The association between a given module and the trait of interest was calculated using bicor values. All modules whose ME displays a significant (FDR adjusted p-values < 0.05), moderate or higher (≥ 0.5) correlation with age were selected for subsequent analyses. Next, for each of the selected modules, module membership (MM) and gene significance (GS) measures were calculated. MM results from correlating the expression of individual genes to the ME, whereas GS corresponds to the absolute value of the correlation between individual genes and the trait of interest. Similar to the previous step, only modules with moderate or higher (≥ 0.5) and significant (p-values < 0.05) correlations were considered to be relevant. Lastly, for each selected module, genes with individual GS > 0.2 and MM > 0.8 were considered to be the most functionally important, i.e. hub genes (as seen in [91–93]). The R package UpSetR (v. 1.4.0)  was used to calculate and visualize the overlap between Trendy, module, and hub genes.
Functional characterization of Trendy-module-hub genes’ overlap
To functionally characterize the gene lists corresponding to the intersection of Trendy, module, and hub genes per tissue, we performed over-representation analysis of GO BPs using the R package clusterProfiler (v. 3.16.0) . Because this package requires NCBI’s Entrez Gene IDs as input, we converted gene symbols into EntrezIDs with the org.Mm.eg.db R package (v. 3.11.4) . We considered each tissue’s expressed genes (after filtering; see Methods – Data pre-processing and normalization) as the universe for over-representation analysis. GO terms with an FDR adjusted p-value less than 0.05 were selected for subsequent analyses.
Network visualization of functionally enriched terms
Network visualization of the enriched GO terms used the enrichmentMap plugin (v. 3.3.0)  of Cytoscape (v. 3.8.0) , with nodes representing GO terms, and edges depicting similarity scores based on the number of genes in common between nodes. To construct our networks, we set an edge similarity cutoff of 0.7. GO term redundancy was addressed with the AutoAnnotate (v. 1.3.3) , clusterMaker2 (v. 1.3.1) , and WordCloud (v. 3.1.3)  plug-ins from Cytoscape. Similar GO terms were clustered together using the Markov Clustering Algorithm (MCL), also with an edge similarity cutoff of 0.7, and cluster labels were created with the default label algorithm Adjacent Words, with 4 maximum words per label and an adjacent word bonus of 8.
All analyses in the R statistical environment were performed using R (v. 4.0.3) and RStudio (v. 1.1.463). All R code is available at https://github.com/mmccferreira/Aging_2021.
AD: Alzheimer’s Disease; ATP: Adenosine triphosphate; Bicor: Biweight Midcorrelation; BP: Biological Process; DEGs: Differentially Expressed Genes; ER: endoplasmic reticulum; FDR: False Discovery Rate; GEO: Gene Expression Omnibus; GO: Gene Ontology; GS: Gene Significance; HCC: hepatocellular carcinoma; KC: Kupffer cells; ME: Module Eigengene; MHCI: Major Histocompatibility Complex I; MM: Module Membership; MS: Multiple Sclerosis; NAFLD: Non-alcoholic fatty liver disease; PCA: Principal Component Analysis; RNA-Seq: RNA Sequencing; Serpina3n: serine (or cysteine) peptidase inhibitor, clade A, member 3N; Serpinf1: serine (or cysteine) peptidase inhibitor, clade F, member 1; Serping1: serine (or cysteine) peptidase inhibitor, clade F, member 1; TOM: Topological Overlap Matrices; VST: Variance Stabilizing Transformation; WGCNA: Weighted Gene Correlation Network Analysis; z.K: Standardized Connectivity.
M.F. carried out the RNASeq data analysis and the gene network analysis, participated in the definition of the questions and wrote the paper; S.F. wrote and corrected the paper; M.P. and A.R. developed pipelines for data analysis; A.N. and A.R.S. corrected the paper; A.J.R., N.S. and G.M. contributed to the definition of the experimental work plan; and M.A.S.S. coordinated the overall project, defined the questions and defined the work plan, wrote and corrected the paper.
We are most thankful to the University of Aveiro Genome Medicine Laboratory and Institute of Biomedicine - iBiMED for supporting this work.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by the Portuguese Foundation for Science and Technology (FCT) and FEDER (Fundo Europeu de Desenvolvimento Regional) funds through the COMPETE 2020, Operational Programme for Competitiveness and Internationalization (POCI) (GenomePT: POCI-01-0145-FEDER-022184; MEDPERSYST: POCI-01-0145-FEDER-016428-PAC; WISDOM: POCI-01-0145-FEDER-029843) and by Centro 2020 program, Portugal 2020 and European Regional Development Fund (pAGE Integrated project: Centro-01-0145-FEDER-000003; MEDISIS: CENTRO-01-0246-FEDER-000018; NETDIAMOND: POCI-01-0145-FEDER-016385). The iBiMED research unit is supported by FCT (UID/BIM/04501/2020). MF and SF are directly supported by FCT grants (SFRH/BD/131736/2017 and SFRH/BD/148323/2019).
- 1. Lee JS, Ward WO, Ren H, Vallanat B, Darlington GJ, Han ES, Laguna JC, DeFord JH, Papaconstantinou J, Selman C, Corton JC. Meta-analysis of gene expression in the mouse liver reveals biomarkers associated with inflammation increased early during aging. Mech Ageing Dev. 2012; 133:467–78. https://doi.org/10.1016/j.mad.2012.05.006 [PubMed]
- 2. Voutetakis K, Chatziioannou A, Gonos ES, Trougakos IP. Comparative Meta-Analysis of Transcriptomics Data during Cellular Senescence and In Vivo Tissue Ageing. Oxid Med Cell Longev. 2015; 2015:732914. https://doi.org/10.1155/2015/732914 [PubMed]
- 3. Benayoun BA, Pollina EA, Singh PP, Mahmoudi S, Harel I, Casey KM, Dulken BW, Kundaje A, Brunet A. Remodeling of epigenome and transcriptome landscapes with aging in mice reveals widespread induction of inflammatory responses. Genome Res. 2019; 29:697–709. https://doi.org/10.1101/gr.240093.118 [PubMed]
- 4. Clark D, Brazina S, Yang F, Hu D, Hsieh CL, Niemi EC, Miclau T, Nakamura MC, Marcucio R. Age-related changes to macrophages are detrimental to fracture healing in mice. Aging Cell. 2020; 19:e13112. https://doi.org/10.1111/acel.13112 [PubMed]
- 5. Zahn JM, Poosala S, Owen AB, Ingram DK, Lustig A, Carter A, Weeraratna AT, Taub DD, Gorospe M, Mazan-Mamczarz K, Lakatta EG, Boheler KR, Xu X, et al. AGEMAP: a gene expression database for aging in mice. PLoS Genet. 2007; 3:e201. https://doi.org/10.1371/journal.pgen.0030201 [PubMed]
- 6. Schaum N, Lehallier B, Hahn O, Pálovics R, Hosseinzadeh S, Lee SE, Sit R, Lee DP, Losada PM, Zardeneta ME, Fehlmann T, Webber JT, McGeever A, et al, and Tabula Muris Consortium. Ageing hallmarks exhibit organ-specific temporal signatures. Nature. 2020; 583:596–602. https://doi.org/10.1038/s41586-020-2499-y [PubMed]
- 7. Glass D, Viñuela A, Davies MN, Ramasamy A, Parts L, Knowles D, Brown AA, Hedman AK, Small KS, Buil A, Grundberg E, Nica AC, Di Meglio P, et al, UK Brain Expression consortium, and MuTHER consortium. Gene expression changes with age in skin, adipose tissue, blood and brain. Genome Biol. 2013; 14:R75. https://doi.org/10.1186/gb-2013-14-7-r75 [PubMed]
- 8. Bakken TE, Miller JA, Ding SL, Sunkin SM, Smith KA, Ng L, Szafer A, Dalley RA, Royall JJ, Lemon T, Shapouri S, Aiona K, Arnold J, et al. A comprehensive transcriptional map of primate brain development. Nature. 2016; 535:367–75. https://doi.org/10.1038/nature18637 [PubMed]
- 9. Palmer D, Fabris F, Doherty A, Freitas AA, de Magalhães JP. Ageing transcriptome meta-analysis reveals similarities and differences between key mammalian tissues. Aging (Albany NY). 2021; 13:3313–41. https://doi.org/10.18632/aging.202648 [PubMed]
- 10. Merienne N, Meunier C, Schneider A, Seguin J, Nair SS, Rocher AB, Le Gras S, Keime C, Faull R, Pellerin L, Chatton JY, Neri C, Merienne K, Déglon N. Cell-Type-Specific Gene Expression Profiling in Adult Mouse Brain Reveals Normal and Disease-State Signatures. Cell Rep. 2019; 26:2477–93.e9. https://doi.org/10.1016/j.celrep.2019.02.003 [PubMed]
- 11. Frenk S, Houseley J. Gene expression hallmarks of cellular ageing. Biogerontology. 2018; 19:547–66. https://doi.org/10.1007/s10522-018-9750-z [PubMed]
- 12. Harries LW, Hernandez D, Henley W, Wood AR, Holly AC, Bradley-Smith RM, Yaghootkar H, Dutta A, Murray A, Frayling TM, Guralnik JM, Bandinelli S, Singleton A, et al. Human aging is characterized by focused changes in gene expression and deregulation of alternative splicing. Aging Cell. 2011; 10:868–78. https://doi.org/10.1111/j.1474-9726.2011.00726.x [PubMed]
- 13. Zhou X, Zhen X, Liu Y, Cui Z, Yue Z, Xu A, Han J. Identification of Key Modules, Hub Genes, and Noncoding RNAs in Chronic Rhinosinusitis with Nasal Polyps by Weighted Gene Coexpression Network Analysis. Biomed Res Int. 2020; 2020:6140728. https://doi.org/10.1155/2020/6140728 [PubMed]
- 14. Cellerino A, Ori A. What have we learned on aging from omics studies? Semin Cell Dev Biol. 2017; 70:177–89. https://doi.org/10.1016/j.semcdb.2017.06.012 [PubMed]
- 15. Gomez-Verjan JC, Vazquez-Martinez ER, Rivero-Segura NA, Medina-Campos RH. The RNA world of human ageing. Hum Genet. 2018; 137:865–79. https://doi.org/10.1007/s00439-018-1955-3 [PubMed]
- 16. de Magalhães JP, Curado J, Church GM. Meta-analysis of age-related gene expression profiles identifies common signatures of aging. Bioinformatics. 2009; 25:875–81. https://doi.org/10.1093/bioinformatics/btp073 [PubMed]
- 17. Stegeman R, Weake VM. Transcriptional Signatures of Aging. J Mol Biol. 2017; 429:2427–37. https://doi.org/10.1016/j.jmb.2017.06.019 [PubMed]
- 18. Anisimova AS, Meerson MB, Gerashchenko MV, Kulakovskiy IV, Dmitriev SE, Gladyshev VN. Multifaceted deregulation of gene expression and protein synthesis with age. Proc Natl Acad Sci USA. 2020; 117:15581–90. https://doi.org/10.1073/pnas.2001788117 [PubMed]
- 19. Aramillo Irizar P, Schäuble S, Esser D, Groth M, Frahm C, Priebe S, Baumgart M, Hartmann N, Marthandan S, Menzel U, Müller J, Schmidt S, Ast V, et al. Transcriptomic alterations during ageing reflect the shift from cancer to degenerative diseases in the elderly. Nat Commun. 2018; 9:327. https://doi.org/10.1038/s41467-017-02395-2 [PubMed]
- 20. Srivastava A, Barth E, Ermolaeva MA, Guenther M, Frahm C, Marz M, Witte OW. Tissue-specific Gene Expression Changes Are Associated with Aging in Mice. Genomics Proteomics Bioinformatics. 2020; 18:430–42. https://doi.org/10.1016/j.gpb.2020.12.001 [PubMed]
- 21. Ori A, Toyama BH, Harris MS, Bock T, Iskar M, Bork P, Ingolia NT, Hetzer MW, Beck M. Integrated Transcriptome and Proteome Analyses Reveal Organ-Specific Proteome Deterioration in Old Rats. Cell Syst. 2015; 1:224–37. https://doi.org/10.1016/j.cels.2015.08.012 [PubMed]
- 22. Zhou Q, Wan Q, Jiang Y, Liu J, Qiang L, Sun L. A Landscape of Murine Long Non-Coding RNAs Reveals the Leading Transcriptome Alterations in Adipose Tissue during Aging. Cell Rep. 2020; 31:107694. https://doi.org/10.1016/j.celrep.2020.107694 [PubMed]
- 23. Fu C, Hickey M, Morrison M, McCarter R, Han ES. Tissue specific and non-specific changes in gene expression by aging and by early stage CR. Mech Ageing Dev. 2006; 127:905–16. https://doi.org/10.1016/j.mad.2006.09.006 [PubMed]
- 24. Melis JP, Jonker MJ, Vijg J, Hoeijmakers JH, Breit TM, van Steeg H. Aging on a different scale--chronological versus pathology-related aging. Aging (Albany NY). 2013; 5:782–88. https://doi.org/10.18632/aging.100606 [PubMed]
- 25. Zhang MJ, Pisco AO, Darmanis S, Zou J. Mouse aging cell atlas analysis reveals global and cell type-specific aging signatures. Elife. 2021; 10:e62293. https://doi.org/10.7554/eLife.62293 [PubMed]
- 26. Almanzar N, Antony J, Baghel AS, Bakerman I, Bansal I, Barres BA, Beachy PA, Berdnik D, Bilen B, Brownfield D, Cain C, Chan CKF, Chen MB, et al. A single-cell transcriptomic atlas characterizes ageing tissues in the mouse. Nat 2020 5837817 [Internet]. Nature Publishing Group; 2020 [cited 2021 Jul 9]; 583:590–5. https://www.nature.com/articles/s41586-020-2496-1.
- 27. Beveridge NJ, Santarelli DM, Wang X, Tooney PA, Webster MJ, Weickert CS, Cairns MJ. Maturation of the human dorsolateral prefrontal cortex coincides with a dynamic shift in microRNA expression. Schizophr Bull. 2014; 40:399–409. https://doi.org/10.1093/schbul/sbs198 [PubMed]
- 28. Colantuoni C, Lipska BK, Ye T, Hyde TM, Tao R, Leek JT, Colantuoni EA, Elkahloun AG, Herman MM, Weinberger DR, Kleinman JE. Temporal dynamics and genetic control of transcription in the human prefrontal cortex. Nature. 2011; 478:519–23. https://doi.org/10.1038/nature10524 [PubMed]
- 29. Erraji-Benchekroun L, Underwood MD, Arango V, Galfalvy H, Pavlidis P, Smyrniotopoulos P, Mann JJ, Sibille E. Molecular aging in human prefrontal cortex is selective and continuous throughout adult life. Biol Psychiatry. 2005; 57:549–58. https://doi.org/10.1016/j.biopsych.2004.10.034 [PubMed]
- 30. Cao K, Ryvkin P, Hwang YC, Johnson FB, Wang LS. Analysis of nonlinear gene expression progression reveals extensive pathway and age-specific transitions in aging human brains. PLoS One. 2013; 8:e74578. https://doi.org/10.1371/journal.pone.0074578 [PubMed]
- 31. Irizar H, Goñi J, Alzualde A, Castillo-Triviño T, Olascoaga J, Lopez de Munain A, Otaegui D. Age gene expression and coexpression progressive signatures in peripheral blood leukocytes. Exp Gerontol. 2015; 72:50–56. https://doi.org/10.1016/j.exger.2015.09.003 [PubMed]
- 32. Timmons JA, Volmar CH, Crossland H, Phillips BE, Sood S, Janczura KJ, Törmäkangas T, Kujala UM, Kraus WE, Atherton PJ, Wahlestedt C. Longevity-related molecular pathways are subject to midlife “switch” in humans. Aging Cell. 2019; 18:e12970. https://doi.org/10.1111/acel.12970 [PubMed]
- 33. Yu Y, Fuscoe JC, Zhao C, Guo C, Jia M, Qing T, Bannon DI, Lancashire L, Bao W, Du T, Luo H, Su Z, Jones WD, et al. A rat RNA-Seq transcriptomic BodyMap across 11 organs and 4 developmental stages. Nat Commun. 2014; 5:3230. https://doi.org/10.1038/ncomms4230 [PubMed]
- 34. When are mice considered old? 2020. https://www.jax.org/news-and-insights/jax-blog/2017/november/when-are-mice-considered-old.
- 35. Jonker MJ, Melis JP, Kuiper RV, van der Hoeven TV, Wackers PF, Robinson J, van der Horst GT, Dollé ME, Vijg J, Breit TM, Hoeijmakers JH, van Steeg H. Life spanning murine gene expression profiles in relation to chronological and pathological aging in multiple organs. Aging Cell. 2013; 12:901–09. https://doi.org/10.1111/acel.12118 [PubMed]
- 36. Yang X, Schadt EE, Wang S, Wang H, Arnold AP, Ingram-Drake L, Drake TA, Lusis AJ. Tissue-specific expression and regulation of sexually dimorphic genes in mice. Genome Res. 2006; 16:995–1004. https://doi.org/10.1101/gr.5217506 [PubMed]
- 37. Lopes-Ramos CM, Chen CY, Kuijjer ML, Paulson JN, Sonawane AR, Fagny M, Platig J, Glass K, Quackenbush J, DeMeo DL. Sex Differences in Gene Expression and Regulatory Networks across 29 Human Tissues. Cell Rep. 2020; 31:107795. https://doi.org/10.1016/j.celrep.2020.107795 [PubMed]
- 38. Naqvi S, Godfrey AK, Hughes JF, Goodheart ML, Mitchell RN, Page DC. Conservation, acquisition, and functional impact of sex-biased gene expression in mammals. Science. 2019; 365:eaaw7317. https://doi.org/10.1126/science.aaw7317 [PubMed]
- 39. Sampathkumar NK, Bravo JI, Chen Y, Danthi PS, Donahue EK, Lai RW, Lu R, Randall LT, Vinson N, Benayoun BA. Widespread sex dimorphism in aging and age-related diseases. Hum Genet. 2020; 139:333–56. https://doi.org/10.1007/s00439-019-02082-w [PubMed]
- 40. Wang Y, Eng DG, Pippin JW, Gharib SA, McClelland A, Gross KW, Shankland SJ. Sex differences in transcriptomic profiles in aged kidney cells of renin lineage. Aging (Albany NY). 2018; 10:606–21. https://doi.org/10.18632/aging.101416 [PubMed]
- 41. Du M, Mangold CA, Bixler GV, Brucklacher RM, Masser DR, Stout MB, Elliott MH, Freeman WM. Retinal gene expression responses to aging are sexually divergent. Mol Vis [Internet]. Molecular Vision; 2017 [cited 2021 Jun 21]; 23: 707–17. http://www.molvis.org/molvis/v23/707.
- 42. Mangold CA, Wronowski B, Du M, Masser DR, Hadad N, Bixler GV, Brucklacher RM, Ford MM, Sonntag WE, Freeman WM. Sexually divergent induction of microglial-associated neuroinflammation with hippocampal aging. J Neuroinflammation. 2017; 14:141. https://doi.org/10.1186/s12974-017-0920-8 [PubMed]
- 43. Oh VS, Li RW. Temporal Dynamic Methods for Bulk RNA-Seq Time Series Data. Genes (Basel). 2021; 12:352. https://doi.org/10.3390/genes12030352 [PubMed]
- 44. GEO Accession viewer. 2020. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132040.
- 45. Tabula Muris Consortium, Overall coordination, Logistical coordination, Organ collection and processing, Library preparation and sequencing, Computational data analysis, Cell type annotation, Writing group, Supplemental text writing group, and Principal investigators. Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature. 2018; 562:367–72. https://doi.org/10.1038/s41586-018-0590-4 [PubMed]
- 46. Bacher R, Leng N, Chu LF, Ni Z, Thomson JA, Kendziorski C, Stewart R. Trendy: segmented regression analysis of expression dynamics in high-throughput ordered profiling experiments. BMC Bioinformatics. 2018; 19:380. https://doi.org/10.1186/s12859-018-2405-x [PubMed]
- 47. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9:559. https://doi.org/10.1186/1471-2105-9-559 [PubMed]
- 48. Barry C, Schmitz MT, Argus C, Bolin JM, Probasco MD, Leng N, Duffin BM, Steill J, Swanson S, McIntosh BE, Stewart R, Kendziorski C, Thomson JA, Bacher R. Automated minute scale RNA-seq of pluripotent stem cell differentiation reveals early divergence of human and mouse gene expression kinetics. PLoS Comput Biol. 2019; 15:e1007543. https://doi.org/10.1371/journal.pcbi.1007543 [PubMed]
- 49. Chen A, Wang Y, Wu J, Tang D, Zhu Q, Lu A, Yang J, Cai Z, Shi J. Identification and characterization of dynamically regulated hepatitis-related genes in a concanavalin A-induced liver injury model. Aging (Albany NY). 2020; 12:23187–99. https://doi.org/10.18632/aging.104089 [PubMed]
- 50. Li C, Mo D, Li M, Zheng Y, Li Q, Ou S, Zhang Z. Age-related but not longevity-related genes are found by weighted gene co-expression network analysis in the peripheral blood cells of humans. Genes Genet Syst. 2019; 93:221–28. https://doi.org/10.1266/ggs.17-00052 [PubMed]
- 51. Hu Y, Pan J, Xin Y, Mi X, Wang J, Gao Q, Luo H. Gene Expression Analysis Reveals Novel Gene Signatures Between Young and Old Adults in Human Prefrontal Cortex. Front Aging Neurosci. 2018; 10:259. https://doi.org/10.3389/fnagi.2018.00259 [PubMed]
- 52. Li ML, Wu SH, Zhang JJ, Tian HY, Shao Y, Wang ZB, Irwin DM, Li JL, Hu XT, Wu DD. Author Correction: 547 transcriptomes from 44 brain areas reveal features of the aging brain in non-human primates. Genome Biol. 2020; 21:46. https://doi.org/10.1186/s13059-020-01962-z [PubMed]
- 53. Ximerakis M, Lipnick SL, Innes BT, Simmons SK, Adiconis X, Dionne D, Mayweather BA, Nguyen L, Niziolek Z, Ozek C, Butty VL, Isserlin R, Buchanan SM, et al. Single-cell transcriptomic profiling of the aging mouse brain. Nat Neurosci. 2019; 22:1696–708. https://doi.org/10.1038/s41593-019-0491-3 [PubMed]
- 54. McKeever PM, Kim T, Hesketh AR, MacNair L, Miletic D, Favrin G, Oliver SG, Zhang Z, St George-Hyslop P, Robertson J. Cholinergic neuron gene expression differences captured by translational profiling in a mouse model of Alzheimer’s disease. Neurobiol Aging. 2017; 57:104–19. https://doi.org/10.1016/j.neurobiolaging.2017.05.014 [PubMed]
- 55. Kucera M, Isserlin R, Arkhangorodsky A, Bader GD. AutoAnnotate: A Cytoscape app for summarizing networks with semantic annotations. F1000Res. 2016; 5:1717. https://doi.org/10.12688/f1000research.9090.1 [PubMed]
- 56. Sayegh JF, Lajtha A. In vivo rates of protein synthesis in brain, muscle, and liver of five vertebrate species. Neurochem Res. 1989; 14:1165–68. https://doi.org/10.1007/BF00965625 [PubMed]
- 57. Shahbazian FM, Jacobs M, Lajtha A. Rates of protein synthesis in brain and other organs. Int J Dev Neurosci. 1987; 5:39–42. https://doi.org/10.1016/0736-5748(87)90046-3 [PubMed]
- 58. Augereau C, Lemaigre FP, Jacquemin P. Extraction of high-quality RNA from pancreatic tissues for gene expression studies. Anal Biochem. 2016; 500:60–62. https://doi.org/10.1016/j.ab.2016.02.008 [PubMed]
- 59. Azevedo-Pouly AC, Elgamal OA, Schmittgen TD. RNA isolation from mouse pancreas: a ribonuclease-rich tissue. J Vis Exp. 2014; e51779. https://doi.org/10.3791/51779 [PubMed]
- 60. InanlooRahatloo K, Liang G, Vo D, Ebert A, Nguyen I, Nguyen PK. Sex-based differences in myocardial gene expression in recently deceased organ donors with no prior cardiovascular disease. PLoS One. 2017; 12:e0183874. https://doi.org/10.1371/journal.pone.0183874 [PubMed]
- 61. Sanchez D, Vallée B, Figarella C, Baeza N, Guy-Crotte O. Sexual dimorphism of pancreatic gene expression in normal mice. Pancreas. 2000; 21:407–13. https://doi.org/10.1097/00006676-200011000-00012 [PubMed]
- 62. Rossetti CL, de Oliveira Costa HM, Barthem CS, da Silva MH, de Carvalho DP, da-Silva WS. Sexual dimorphism of liver endoplasmic reticulum stress susceptibility in prepubertal rats and the effect of sex steroid supplementation. Exp Physiol. 2019; 104:677–90. https://doi.org/10.1113/EP087518 [PubMed]
- 63. Jenkins EC, Shah N, Gomez M, Casalena G, Zhao D, Kenny TC, Guariglia SR, Manfredi G, Germain D. Proteasome mapping reveals sexual dimorphism in tissue-specific sensitivity to protein aggregations. EMBO Rep. 2020; 21:e48978. https://doi.org/10.15252/embr.201948978 [PubMed]
- 64. Francisco S, Ferreira M, Moura G, Soares AR, Santos MA. Does proteostasis get lost in translation? Implications for protein aggregation across the lifespan. Ageing Res Rev. 2020; 62:101119. https://doi.org/10.1016/j.arr.2020.101119 [PubMed]
- 65. López-Otín C, Blasco MA, Partridge L, Serrano M, Kroemer G. The hallmarks of aging. Cell. 2013; 153:1194–217. https://doi.org/10.1016/j.cell.2013.05.039 [PubMed]
- 66. Walther DM, Kasturi P, Zheng M, Pinkert S, Vecchi G, Ciryam P, Morimoto RI, Dobson CM, Vendruscolo M, Mann M, Hartl FU. Widespread Proteome Remodeling and Aggregation in Aging C. elegans. Cell. 2015; 161:919–32. https://doi.org/10.1016/j.cell.2015.03.032 [PubMed]
- 67. David DC, Ollikainen N, Trinidad JC, Cary MP, Burlingame AL, Kenyon C. Widespread protein aggregation as an inherent part of aging in C. elegans. PLoS Biol. 2010; 8:e1000450. https://doi.org/10.1371/journal.pbio.1000450 [PubMed]
- 68. Hamer G, Matilainen O, Holmberg CI. A photoconvertible reporter of the ubiquitin-proteasome system in vivo. Nat Methods. 2010; 7:473–78. https://doi.org/10.1038/nmeth.1460 [PubMed]
- 69. Groh N, Bühler A, Huang C, Li KW, van Nierop P, Smit AB, Fändrich M, Baumann F, David DC. Age-Dependent Protein Aggregation Initiates Amyloid-β Aggregation. Front Aging Neurosci. 2017; 9:138. https://doi.org/10.3389/fnagi.2017.00138 [PubMed]
- 70. Chen X, Miller NM, Afghah Z, Geiger JD. Development of AD-Like Pathology in Skeletal Muscle. J Parkinsons Dis Alzheimers Dis. 2019; 6. https://doi.org/10.13188/2376-922x.1000028 [PubMed]
- 71. Askanas V, Engel WK, Nogalska A. Sporadic inclusion-body myositis: A degenerative muscle disease associated with aging, impaired muscle protein homeostasis and abnormal mitophagy. Biochim Biophys Acta. 2015; 1852:633–43. https://doi.org/10.1016/j.bbadis.2014.09.005 [PubMed]
- 72. Babagana M, Oh KS, Chakraborty S, Pacholewska A, Aqdas M, Sung MH. Hedgehog dysregulation contributes to tissue-specific inflammaging of resident macrophages. bioRxiv [Internet]. 2020. https://doi.org/10.1101/2020.08.16.253237
- 73. Lucas A, Yaron JR, Zhang L, Ambadapadi S. Overview of Serpins and Their Roles in Biological Systems. Methods Mol Biol. 2018; 1826:1–7. https://doi.org/10.1007/978-1-4939-8645-3_1 [PubMed]
- 74. Ho TC, Chiang YP, Chuang CK, Chen SL, Hsieh JW, Lan YW, Tsao YP. PEDF-derived peptide promotes skeletal muscle regeneration through its mitogenic effect on muscle progenitor cells. Am J Physiol Cell Physiol. 2015; 309:C159–68. https://doi.org/10.1152/ajpcell.00344.2014 [PubMed]
- 75. Okada K, Naito AT, Higo T, Nakagawa A, Shibamoto M, Sakai T, Hashimoto A, Kuramoto Y, Sumida T, Nomura S, Ito M, Yamaguchi T, Oka T, et al. Wnt/β-Catenin Signaling Contributes to Skeletal Myopathy in Heart Failure via Direct Interaction With Forkhead Box O. Circ Heart Fail. 2015; 8:799–808. https://doi.org/10.1161/CIRCHEARTFAILURE.114.001958 [PubMed]
- 76. Kamei Y, Miura S, Suzuki M, Kai Y, Mizukami J, Taniguchi T, Mochida K, Hata T, Matsuda J, Aburatani H, Nishino I, Ezaki O. Skeletal muscle FOXO1 (FKHR) transgenic mice have less skeletal muscle mass, down-regulated Type I (slow twitch/red muscle) fiber genes, and impaired glycemic control. J Biol Chem. 2004; 279:41114–23. https://doi.org/10.1074/jbc.M400674200 [PubMed]
- 77. Boisvert MM, Erikson GA, Shokhirev MN, Allen NJ. The Aging Astrocyte Transcriptome from Multiple Regions of the Mouse Brain. Cell Rep. 2018; 22:269–85. https://doi.org/10.1016/j.celrep.2017.12.039 [PubMed]
- 78. Vanni S, Moda F, Zattoni M, Bistaffa E, De Cecco E, Rossi M, Giaccone G, Tagliavini F, Haïk S, Deslys JP, Zanusso G, Ironside JW, Ferrer I, et al. Differential overexpression of SERPINA3 in human prion diseases. Sci Rep. 2017; 7:15637. https://doi.org/10.1038/s41598-017-15778-8 [PubMed]
- 79. Aslam MS, Yuan L. Serpina3n: Potential drug and challenges, mini review. J Drug Target. 2020; 28:368–78. https://doi.org/10.1080/1061186X.2019.1693576 [PubMed]
- 80. Fissolo N, Matute-Blanch C, Osman M, Costa C, Pinteac R, Miró B, Sanchez A, Brito V, Dujmovic I, Voortman M, Khalil M, Borràs E, Sabidó E, et al. CSF SERPINA3 Levels Are Elevated in Patients With Progressive MS. Neurol Neuroimmunol Neuroinflamm. 2021; 8:e941. https://doi.org/10.1212/NXI.0000000000000941 [PubMed]
- 81. Tabula Muris Consortium, Hosseinzadeh S. SmartSeq2 for HTP Generation of Bulk RNA Libraries V.1. 2019. https://doi.org/10.17504/protocols.io.2uvgew6
- 82. Oldham MC, Langfelder P, Horvath S. Network methods for describing sample relationships in genomic datasets: application to Huntington’s disease. BMC Syst Biol. 2012; 6:63. https://doi.org/10.1186/1752-0509-6-63 [PubMed]
- 83. Horvath S. Corrected R code from chapter 12 of the book 2020. http://pages.stat.wisc.edu/~yandell/statgen/ucla/WGCNA/wgcna.html.
- 84. Yates AD, Achuthan P, Akanni W, Allen J, Allen J, Alvarez-Jarreta J, Amode MR, Armean IM, Azov AG, Bennett R, Bhai J, Billis K, Boddu S, et al. Ensembl 2020. Nucleic Acids Res. 2020; 48:D682–88. https://doi.org/10.1093/nar/gkz966 [PubMed]
- 85. Durinck S, Moreau Y, Kasprzyk A, Davis S, De Moor B, Brazma A, Huber W. BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics. 2005; 21:3439–40. https://doi.org/10.1093/bioinformatics/bti525 [PubMed]
- 86. Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc. 2009; 4:1184–91. https://doi.org/10.1038/nprot.2009.97 [PubMed]
- 87. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550. https://doi.org/10.1186/s13059-014-0550-8 [PubMed]
- 88. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11:R106. https://doi.org/10.1186/gb-2010-11-10-r106 [PubMed]
- 89. Tutorials for WGCNA R package. 2020. https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/index.html.
- 90. WGCNA package: Frequently Asked Questions. 2020. https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html.
- 91. Liu Y, Gu HY, Zhu J, Niu YM, Zhang C, Guo GL. Identification of Hub Genes and Key Pathways Associated With Bipolar Disorder Based on Weighted Gene Co-expression Network Analysis. Front Physiol. 2019; 10:1081. https://doi.org/10.3389/fphys.2019.01081 [PubMed]
- 92. Wang W, Jiang W, Hou L, Duan H, Wu Y, Xu C, Tan Q, Li S, Zhang D. Weighted gene co-expression network analysis of expression data of monozygotic twins identifies specific modules and hub genes related to BMI. BMC Genomics. 2017; 18:872. https://doi.org/10.1186/s12864-017-4257-6 [PubMed]
- 93. Yuan L, Chen L, Qian K, Qian G, Wu CL, Wang X, Xiao Y. Co-expression network analysis identified six hub genes in association with progression and prognosis in human clear cell renal cell carcinoma (ccRCC). Genom Data. 2017; 14:132–40. https://doi.org/10.1016/j.gdata.2017.10.006 [PubMed]
- 94. Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017; 33:2938–40. https://doi.org/10.1093/bioinformatics/btx364 [PubMed]
- 95. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 96. Carlson M. org.Mm.eg.db: Genome wide annotation for Mouse. R package version 3.8.2. Bioconductor. 2019. http://bioconductor.org/packages/release/data/annotation/html/org.Mm.eg.db.html.
- 97. Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS One. 2010; 5:e13984. https://doi.org/10.1371/journal.pone.0013984 [PubMed]
- 98. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
- 99. Morris JH, Apeltsin L, Newman AM, Baumbach J, Wittkop T, Su G, Bader GD, Ferrin TE. clusterMaker: a multi-algorithm clustering plugin for Cytoscape. BMC Bioinformatics. 2011; 12:436. https://doi.org/10.1186/1471-2105-12-436 [PubMed]
- 100. Oesper L, Merico D, Isserlin R, Bader GD. WordCloud: a Cytoscape plugin to create a visual semantic summary of networks. Source Code Biol Med. 2011; 6:7. https://doi.org/10.1186/1751-0473-6-7 [PubMed]