A combination of differential expression and network connectivity analyses identifies a common set of RNA splicing and processing genes altered with age across human tissues

Dataset

Gene expression data (in TPM – transcripts per million) were obtained from the Genotype-Tissue Expression (GTEx) project. GTEx (v8 release) has postmortem RNA-seq data from 948 donors aged 20 to 70, collected from 54 tissues, totaling 17,382 samples. To increase the power of the analyses and allow the division of the dataset into age ranges, with an adequate number of samples in each age range, we selected only tissues with at least 800 samples (blood, brain, adipose tissue, muscle, blood vessel, heart, skin, and esophagus). Tissues were analyzed at the level defined by the SMTS sample attribute in the GTEx annotation files.

Differential expression analysis with linear mixed model

Let Y_ijk be the expression level (log2 normalized TPM + 1) of gene j in sample i of individual k, Age_i be the age of sample i, Sex_i be the sex of sample i, Race_i be the race of sample i, Ischemia_i be the time of ischemia of sample i, Batch_i be the experimental batch in which sample i’s gene expression was measured, Death_i be the type of death, in the Hardy scale, of the donor from which sample i was collected, and ε_ijk be the error term (assumed to be normally and independently distributed with var(ε) = σ²I). Then, we identified differentially expressed genes for each tissue separately using the following linear mixed model:

$$\begin{array}{l}{Y}_{ijk}={\beta }_{0j}+{\beta }_{1j}Ag{e}_i+{\beta }_{2j}Se{x}_i+{\beta }_{3j}Rac{e}_i\\ +{\beta }_{4j}Ischemi{a}_i+{\beta }_{5j}Batc{h}_i\\ +{\beta }_{6j}Deat{h}_i+{\beta }_{7j}Individua{l}_k\\ +{\in }_{ijk}\end{array}$$

Individuals are considered a block random effect, considering that some tissues have more than one sample for the same individual. Ischemic time was divided into 300-minute intervals to standardize the procedure with the removal of confounding factors, which was done in a later step. Random effects (Individual_k) are assumed to be normally and independently distributed with var(Individual) = σ_ind²I.

Genes with no or slight variations in expression were filtered out of the analysis. Selected genes had expression higher than 0.1 TPM in at least 20% of the samples. Depending on the analysis, genes needed to satisfy an FDR (false discovery rate) [19] <0.1 or <0.25 to be considered differentially expressed genes (DEGs) with age, considering BH multiple test corrections. The higher threshold of 0.25 was considered due to the rationale detailed in Results - Comparisons of altered genes across tissues. The significance of these higher threshold results is validated by the permutation tests described in Methods - Permutation tests for the intersection sets. The model was fit using the nlme package in R [20].

Removal of confounding factors

For network analyses, expression was corrected for confounding factors using ComBat from the sva R package [21]. Since ComBat can only correct for one effect at a time, effects were removed iteratively, following the order of sex, race, experimental batch, ischemic time, and death type. Since ComBat does not accept continuous variables, ischemic time was divided into 300-minute intervals.

Independence of connectivity and mean expression levels

The connectivity of a gene in a co-expression network is defined as the sum of that gene’s co-expression values with every other gene in the network. In co-expression networks, the co-expression between a pair of genes is defined as the correlation in their gene expression levels across samples (elevated to a soft-thresholding power in the case of WGCNA).

Let x_i and y_i be the expression levels of genes X and Y in sample i; $\overline{x}$ and $\overline{y}$ be the mean expression levels of genes X and Y across samples, and n be the total number of samples. Then,

$$cov(X,Y)=\frac{{\displaystyle {\sum }_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}}{n-1}$$

is the covariance between X and Y, and

$${\sigma }_X=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(x_i-\overline{x})}^2}}{n-1}}\text{and}{\sigma }_Y=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{(y_i-\overline{y})}^2}}{n-1}}$$

are the standard deviations of X and Y. Finally, Pearson’s correlation coefficient (r) between a pair of gene expression vectors X and Y is defined as:

$$r_{X,Y}=\frac{cov(X,Y)}{{\sigma }_X{\sigma }_Y}$$

Consider a case in which gene Y is differentially expressed between two conditions, and a value of c increases its mean. If the increase of c is observed uniformly across all samples, the components in the formulas related to Y suffer the following transformation:

$$y_i-\overline{y}\stackrel{+c}{\to }(y_i+c)-(\overline{y}+c)=y_i-\overline{y}$$

Therefore, if the increase in expression is uniform across all samples, the covariances and standard deviations are unaffected. This way, the correlation between X and Y can remain the same in two conditions, even if Y is differentially expressed between the two conditions.

Similarly, correlations between X and Y can change, even if the mean expression values remain constant. Consider two extreme cases:

1. Y = X, with $\overline{y}=\overline{x}=\text{c}:$ all values of y in each sample i are identical to the values of x in the same samples, and means equal c;

2. Y = c, $\overline{y}=\text{c}:$ all Y values in each sample are identical and, consequently, equal to their mean c.

In case 1, the values of Y in each sample are variable and have a mean of c. Since Y and X are identical, their correlation is one. In case 2, Y values are constant and equal to c in each sample. Since Y remains constant as X changes, their correlation is zero. Less extreme cases can be easily imagined.

Therefore, gene expression correlations between pairs of genes can vary significantly between conditions, even if their mean expression values remain unchanged.

Differential connectivity analysis

The dataset was divided into nine 6-year age ranges, except for the last one, which had a 3-year range. This division ensures that the smallest sub-dataset has at least 20 samples, the minimum recommended for WGCNA. Connectivity and eigenvector centrality for each gene were then calculated in each age range using the igraph R package [22]. For comparison between networks, the connectivities inside each age range were ranked, and a linear regression was applied to detect differences in ranked connectivities along the nine age ranges in the two metrics. Depending on the analysis, genes were considered differentially connected genes (DCGs) or differentially eigenconnected genes (DECGs) with FDR<0.1 or FDR<0.25. The higher threshold of 0.25 was considered due to the rationale detailed in Results - Comparisons of altered genes across tissues. The significance of these higher threshold results is validated by the permutation tests described in Methods - Permutation tests for the intersection sets.

Permutation tests for the intersection sets

Intersection sets were tested for significance using permutation tests.

Consider G_ijk as the number of significantly altered genes in tissue i according to regression j (differential expression, differential connectivity, or differential eigenconnectivity) using an FDR threshold k (0.1 or 0.25). We constructed 100,000 permuted sets p (without replacement) of genes with the same size as each G_ijk.

For example, for brain, we sampled without replacement 100,000 sets of 530 genes (number of DEGs at FDR<0.1), 100,000 sets of 7,699 genes (number of DCGs at FDR<0.1), and 100,000 sets of 9,421 genes (number of DECGs at FDR<0.1); with additional three sets for the cases with FDR<0.25.

The same process was repeated for all eight tissues, generating 100,000 sets of genes for each tissue, each of the three regressions, and each FDR threshold (8 tissues x 3 regressions x 2 FDR thresholds x 100,000 permutations). Therefore, each permutation p generated 48 sets, equivalent to all G_ijk.

For each p, we calculated the intersections between tissues, considering all the scenarios analyzed. The proportion of permutations p with intersecting genes higher than the observed values was used to calculate a permutation p-value.

The higher threshold of 0.25 was considered due to the rationale detailed in Results - Comparisons of altered genes across tissues. However, despite this higher FDR threshold chosen for filtering tissue-specific altered genes, all intersection sets chosen for downstream analyses were significant at FDR < 0.001 in the permutation tests.

Co-expression network construction and consensus module identification with weighted gene co-expression network analysis

Gene expression data for each analyzed tissue and age range without the effects of known confounding factors were used as inputs for the WGCNA package in R [23]. WGCNA constructs weighted gene co-expression networks using a soft threshold and fitting the data into a scale-free model. For each tissue, we selected the soft-thresholding power that resulted in networks in all age ranges within that tissue achieving an approximate scale-free topology index of 0.8. The chosen soft-thresholding powers for each tissue were: 6 (blood), 6 (brain), 9 (adipose tissue), 12 (muscle), 8 (blood vessel), 12 (heart), 6 (skin), and 6 (esophagus). We used the consensusBlockwiseModules function for module detection, with networkType set as “unsigned”. All parameters were kept at their default values. Sub-modules were defined by overlapping the positive and negative portions of the sets of genes of the three metrics (positive and negative DEGs, DCGs, and DECGs) with the modules defined by consensus modules analysis.

Functional enrichment analysis

All enrichment analyses were carried out with the topGO R package [24], using the classic Fisher algorithm. Annotations for each gene were retrieved from Ensemble 110 [25]. Multiple test correction was performed using the BH method. Terms were considered enriched with FDR < 0.1.

For the enrichment of intersection sets, when the number of terms was too large to display correctly in the chart, we filtered by only showing the Biological Process terms and, if necessary, removing the terms with fewer genes annotated.

Cross-tissue analysis of connectivity changes

To analyze how the same genes behave in different tissues regarding their connectivity, we only included the 14,489 protein-coding genes present in the analyses of the eight tissues (i.e., those that were not filtered out in any tissue during the initial data filtering stage). To facilitate cross-tissue comparison, we first reranked connectivities using only the protein-coding genes employed in this analysis (i.e., we removed all other genes and recalculated a new ranking in each tissue, including only the 14,489 genes analyzed here).

Our initial ranking included non-protein-coding genes in the calculation to account for their contribution to the information. Common, unfiltered non-protein-coding genes across tissues could have been included in this new ranking. However, given that we are only interested in protein-coding genes, we excluded them from the new rankings because non-protein-coding genes represent around half of the genes available in the dataset, and their connectivities are generally small compared to those of protein-coding genes, being heavily concentrated in the first half of the connectivity rankings. This feature caused significant distortions in the histograms, making visualization and interpretation of the behavior of protein-coding genes (our primary interest) difficult.

We created reference histograms for each tissue, which are histograms of the connectivity ranks of the DCGs for that tissue (referred to as the reference tissue). Given an x-axis representing all the connectivity ranks in bins of 100 ranks, the reference histograms show the connectivity distribution of all genes analyzed. Four reference histograms were constructed for each tissue: positive and negative DCGs in age range one and positive and negative DCGs in age range nine.

After constructing the reference histograms, other histograms were constructed. For each gene considered a DCG in a reference tissue, the connectivity ranks of these genes in different tissues (referred to as target tissues) are presented, regardless of whether they are also DCGs in these other tissues. In cases where DCGs are present in both the reference tissue and the target tissue, the regions representing them in the histograms are colored (red for positive DCGs and blue for negative DCGs).

Changes in correlations were tested for significance using the cocor R package [26], with the cocor.indep.groups function, using Fisher’s z.

Comparison of enriched terms across sub-modules

To compare the results of the enrichment analyses between different sub-modules, we listed all enriched terms with FDR < 0.1 in any sub-module. We constructed a binary matrix with enriched terms in the rows and sub-modules in the columns. Entries had a value of 1 in the matrix if a term was enriched in the respective submodule. The terms were then filtered. First, terms found in fewer than four tissues were removed to focus on processes common across tissues, resulting in a final list of 2,310 terms. Secondly, terms found in the first four and last six levels of the ontology tree from the Biological Process Gene Ontology tree were removed, with the rationale that this strategy removed terms that were too general or too specific to yield relevant biological information. Finally, we removed all the columns (sub-modules) that displayed no entries after the previous filtering. After filtering, 1,956 terms were included in the matrix. These steps ensure that: i) our matrix of terms per sub-module is less sparse (with a higher density of 1 values), making it easier to cluster, interpret, and present (due to its smaller size); and ii) reduce the number of terms to a smaller value that is within the limit of the tool we used for the analysis of semantic similarity (Revigo [27]), which accepts at most 2,000 terms as input.

This binary matrix was then clustered into 8 clusters using the hclust function from the stats package in R [28], using the Ward.D clustering algorithm [29]. The results were presented as a heatmap using the heatmap.2 function from the gplots package in R. The same list of enriched terms that comprise the heatmap rows was then summarized using the online tool Revigo [27]. Revigo summarizes lists of Gene Ontology terms based on their semantic similarities and creates a chart where the distances between points represent their semantic distances. Next to the terms, we included the number of tissues each term is enriched in (considering any sub-module in that tissue). Our list of terms was summarized using the following parameters in Revigo: SimRel semantic similarity algorithm, “Small” option for the size of the final list, “Yes” for removing obsolete terms, “Homo sapiens” for the species, and “Higher value is better” for the interpretation of the additional value (in our case, the number of tissues in which the enriched term is present). The R code for reproducing the chart was then downloaded, and, within R, a variable containing the ID of the cluster in which the hclust function clustered that term was added. The final chart was plotted using the number of tissues as a variable to define the size of the points and the cluster variable to define the coloring of the points.

Protein-protein interaction (PPI) network analysis with STRING

The PPI networks were analyzed using STRING [30], which utilizes the STRING extension for Cytoscape [31]. For constructing our networks, we used scores of at least 0.7, which are considered high in STRING. We queried genes in the 427 set and the 134 set not already included in the 427 set (these sets represent the sets of genes that change in common across 7 and 8 tissues, and how they were selected is described in Results – Comparison of altered genes across tissues), for a total of 506 genes, of which 505 were identified in the STRING database.

Gene overlaps

Gene overlaps and overlap significance statistics were calculated using the GeneOverlap R package [32].

Data availability

The raw data used in this study’s analyses are available through the Genotype-Tissue Expression (GTEx) project at gtexportal.org/home.

Detection of tissue-specific differentially expressed and connected genes

For this study, we used data from the GTEx project (version 8), which has whole transcriptome data for up to 54 tissues from 948 postmortem donors. To maximize statistical power, we aggregated sub-tissues and selected only tissues with at least 800 samples each. We then analyzed the blood, brain, adipose tissue, muscle, blood vessels, heart, skin, and esophagus data. The data was processed as described in Methods. Figure 1A shows a general schematic of data processing and analysis.

The differential expression and connectivity analyses were then carried out as described in the Methods section. For the differential connectivity analyses, the dataset was divided into nine age ranges (Supplementary Figure 1), and networks were constructed using WGCNA [23] (Methods). All results obtained from the three regressions are presented in Supplementary Data 1. For all metrics, genes were considered altered with FDR < 0.1.

The number of differentially expressed genes (DEGs), differentially connected genes (DCGs), and differentially eigenconnected genes (DECGs) varied considerably by tissue (Figure 1B). Notably, many tissues demonstrated relative stability in either the differential expression analysis or one of the connectivity analyses. Blood, brain, muscle, and blood vessels all clearly fall into this category, having many DEGs or DCGs/DECGs, but not both simultaneously. Adipose tissue had many significant genes in all analyses, but a clear tendency towards having more DEGs than DCGs can be seen. Only the heart, skin, and esophagus had relatively similar numbers of DEGs and DCGs/DECGs (though very few in the heart).

Only a small proportion of genes were simultaneously DEGs and DCGs/DECGs (Figure 1C). The p-value distribution of the expression and connectivity regressions against one another shows a significant discordance between the metrics (Supplementary Figure 2A). A similar behavior is seen when comparing the expression regression with the eigenconnectivity regression (Supplementary Figure 2B). On the other hand, the connectivity and eigenconnectivity regressions show remarkable concordance (Supplementary Figure 2B), suggesting these metrics measure similar biological effects.

These results demonstrate that even in tissues with a high number of significant genes across all metrics—such as adipose tissue, skin, and the esophagus—there is a clear tendency for the genes that alter gene expression levels to remain unchanged in one of the connectivity metrics, and vice versa.

Comparisons of altered genes across tissues

To compare the genes found in each tissue, we counted the number of tissues in which each gene was changed in each metric. We define sets of genes changing in common across different tissues as intersection sets. The results for DEGs and DCGs are shown in Figure 2A.

Most intersection sets are significant considering permutation tests (see Methods). Non-significant sets are either in very few tissues (DEGs in 1 and 2 tissues and DCGs in 1 tissue) or in tiny sets in many tissues (DCGs in 5 and 6 tissues). Results, including DECGs, are similar (Supplementary Figure 3).

The non-significance for the sets in one tissue is expected, as they are not true intersections and only indicate that most genes are altered in at least one tissue, demonstrating the heterogeneous component of tissue aging.

Given the discordant nature of the differential expression and differential connectivity regressions (Supplementary Figure 2A), we asked whether genes that are DEGs in some tissues might be DCGs or DECGs in other tissues. This is relevant because both analyses detect changes in transcriptomic behavior, even if the changes are different and independent (see Methods section for a demonstration of the independence between the two metrics). We then considered four scenarios:

1. Whether the gene changed in any metric (DEG/DCG/DECG);

2. Whether the gene is a DEG and/or a DCG;

3. Whether the gene is a DCG and/or a DECG or

4. Whether the gene is a DEG and/or a DECG.

The results show that although most genes changing with age are tissue-specific or shared between a few tissues, when intersections between tissues are considered using different metrics to assess transcriptomic behavior, a core of common genes changes across tissues (Figures 2A and Supplementary Figure 3).

Interestingly, gene sets in scenarios including DEGs and one of the connectivity metrics (DCGs or DECGs) in many tissues were significant, including in eight tissues. Sets that include DCGs and DECGs simultaneously, however, were insignificant. This can be understood considering that DCGs and DECGs overlapped significantly (Supplementary Figure 2B), which makes considering both partially redundant, while permuted sets would be more significant (see Methods).

The small intersection set in eight tissues is caused by the heart. Interestingly, the distribution of raw p-values in the heart does not match what would be expected if the null hypothesis were true for all genes (i.e., none are changing) (Supplementary Figure 4), as can be observed for blood in the differential expression regression (Supplementary Figure 4).

While a typical scenario of FDR correction with so many tests when there are close to no significant genes can be observed in the FDR density histogram for blood in the differential expression regression, the same density histogram in the heart displays a very different pattern, with most FDRs clustering around 0.20-0.25. Additionally, when we plot the FDR distribution of genes belonging to the three sets encompassing more tissues using FDR < 0.1, they cluster around the 0.20-0.25 range in the heart (Supplementary Figure 5). This pattern strongly suggests a non-random distribution of p-values and raises the question of whether the genes that cluster around 0.20-0.25 FDR in the heart are meaningful.

To analyze this, we considered two additional scenarios for the intersection analysis, where we raised the FDR threshold to 0.25:

(i)Considering genes significant at FDR < 0.25 only for the heart and

(ii)Considering FDR < 0.25 for all tissues.

It is essential to note that while increasing the FDR threshold increases the probability of including tissue-specific false positives, this shortcoming is mitigated in the next step of the analysis, where we assess the significance of the intersection sets. Essentially, an actual false positive would be a completely random gene that was below the FDR threshold due to chance alone. However, if the same gene was below the FDR threshold (0.10 or 0.25) in 7 or 8 tissues, then it is extremely unlikely that we are observing a false positive. In other words, true randomness would cause random genes to cross the threshold in different tissues, resulting in very low cross-tissue intersections. So, if a sufficiently large number of genes are crossing the threshold in 7 or 8 tissues, then they are extremely unlikely to be random false positives.

To put it in a slightly different way, any gene that was a false positive in one tissue due to a higher FDR threshold would likely not end up in the final intersection sets, especially in those intersection sets in 7 or 8 tissues, since it would have to coincidentally also be a false positive in all of these other tissues. In essence, the intersection sets in a high number of tissues ultimately exclude false positives. This means that our goal here should be to test whether the number of genes in the intersection sets is large enough to rule out the probability of them being coincidentally false positives.

Indeed, both additional scenarios yielded significant sets in 5, 6, 7, and 8 tissues (Figure 2A); of particular interest are the sets in 8 tissues for both scenarios, with 49 and 427 genes, which represent conserved sets of genes that display transcriptomic alterations across all tested tissues during aging. Both were significant at FDR < 0.001 in the permutation tests, indicating that despite using a higher FDR threshold (0.25) to filter tissue-specific results, the resulting cross-tissue sets are highly significant.

Since the results for the connectivity and eigenconnectivity analyses were well correlated, we focus on DEGs and DCGs in most analyses. Henceforth, we will refer to the sets of genes presented in Figure 2A as “intersection sets,” the specific sets of interest (with the combination of DEGs and DCGs) will be called 3-set, 49-set, 427-set, and 134-set. Figure 2B shows a Venn diagram overlapping these sets of interest.

GO enrichment analysis of intersection sets

We performed a GO enrichment analysis of the intersection sets using the topGO R package. Given that we are not interested in tissue-specific changes but rather in intersection sets across most tissues, we limit the results to intersection sets comprising fewer than 5,000 genes, as shown in Figure 2A. Supplementary Data 2–4 contain enrichment results for all these sets.

In Figure 3, we show the results for selected intersection sets. The sets shown are those with the most tissues for which enrichment results identified significant terms in the Biological Process ontology, accompanied by a reasonable number of essential genes annotated.

While the DEGs chart displays immune system-related terms not found in the DCGs chart, they both contain terms related to RNA splicing and DNA damage response. Furthermore, while the set of DCGs in 4 tissues was not chosen because of the previously mentioned criteria, its analysis shows enrichment in Cellular Component terms for the spliceosomal complex (Supplementary Data 2). This indicates that four tissues display good evidence for enriching intersecting DEGs and DCGs related to RNA splicing and DNA damage response.

The selected intersection sets in seven or eight tissues show similar results, with enrichment in RNA splicing, RNA processing, and response to DNA damage. The FDR distribution of all RNA splicing and processing genes corroborates the tendency of these genes to be located in the areas of greatest significance in the charts of all tissues. It reinforces the tissue-specific discordant behavior of the genes in the two regressions (Supplementary Figure 6).

It is important to note that the term “RNA splicing” is a direct child of the term “RNA processing” in the Biological Process ontology tree, meaning that RNA splicing genes are contained in the RNA processing set.

DNA damage response is not found in the 49-set but is in the 427-set, along with terms related to protein catabolism, autophagy, and ribosome biogenesis. Supplementary Figure 7 displays the induced GO subgraph of the top 20 enriched terms in the 134-set.

In most tissues, genes in these sets are predominantly upregulated in at least one of these two metrics (Supplementary Figure 8), except for blood vessels, which have a mixed profile, and the heart, in which most genes are negative DEGs.

Furthermore, negative DEGs, including those in the 427-set, tended to be highly expressed in younger samples (Supplementary Figures 9, 10). In fact, in blood vessels and the heart, where there was a higher proportion of negative DEGs, both overall and in the intersection sets, the negative DEGs remained highly expressed in older samples, even with the decreases, indicating that most negative DEGs remain necessary for cellular function.

Changes in connectivity between young and old samples tended to be more dramatic, with DCGs in the 427-set and the splicing genes tending to become more highly connected at older ages (Supplementary Figures 9, 10). The most apparent exception seems to be skin, where most splicing and 427-set genes decrease their connectivities (Supplementary Figures 8, 9). The 134-set and the 49-set show similar results (Supplementary Figure 10).

Cross-tissue analysis of connectivity changes

It is puzzling that some tissues tend to have more DEGs or DCGs but not both. Even in the skin and esophagus, where the number of DEGs and DCGs was more balanced, the genes changing in one metric did not change in the other. To explore this phenomenon, we analyzed the distribution of DCG connectivity ranks in each tissue against the connectivity ranks of the same genes in the other tissues (see Methods).

The results indicate that positive DCGs in each tissue tend to be among the lowest or middle-ranked genes in their tissues in younger samples (Figure 4A). However, when we examine the connectivity rank histogram of these same genes in other tissues, we observe a contrasting scenario, with the same set of genes being significantly more highly ranked in the different tissues. For example, positive DCGs in the blood are among the low and middle ranks regarding connectivity. However, the same genes that are positive DCGs in blood tend to be more highly ranked in other tissues.

Negative DCGs show similar behavior with inverse logic. The only exception is skin, since its negative DCGs are highly ranked in many other tissues.

In older samples, the DCGs changed their rankings as expected, such that the connectivity ranking pattern of these genes more closely resembles the pattern they already displayed in other tissues. This suggests that gene connectivity is reorganized in some tissues, and that in tissues where these genes are not changing their connectivities, it is likely because their connectivities were already roughly at the physiologically ideal level.

This convergence of connectivity ranks is evident when comparing the correlations of connectivity ranks between pairs of tissues in the younger and older age ranges. Correlations are consistently significantly increased in all tissues, except correlations involving the skin, which are consistently considerably decreased (Figure 4B – first panel).

A similar pattern can be seen when considering only the intersection sets or all RNA processing and RNA splicing genes (Figure 4B—second to fifth panels). In these cases, however, the reduced number of genes makes it more challenging to determine the significance of all comparisons, especially for the 49-set, even when the increase in correlation is significant (for example, see the comparison between brain and muscle in the 49-set).

Consensus network and GO enrichment analysis of age-related sub-modules

While only a few common genes appear to be changing across 7 or 8 tissues, we inquired whether other genes changing in fewer tissues could exhibit similarities at the process level.

We used WGCNA to build an unsigned consensus network for each tissue’s nine age ranges. This allowed us to identify conserved gene co-expression modules that were present in all nine datasets. After obtaining these modules, we overlapped them with DEGs, DCGs, and DECGs and divided them into positive and negative portions. This led to the definition of sub-modules, where each sub-module comprises the positive or negative DEGs, DCGs, or DECGs contained in that module.

By dividing the DEGs, DCGs, and DECGs using the modular structure provided by WCGNA, we group them into sets likely to be related to the same biological processes and, therefore, identify them based on a functional enrichment analysis.

Since enriching these sub-modules generated many lists of terms, we first summarized them in a heatmap format (Figure 5A and Supplementary Figure 11). Then, we used Revigo to display the results (Figure 5B) (see Methods for a detailed description). Complete enrichment results per module can be found in Supplementary Data 5–7.

The green rows/terms cluster contains sub-modules enriched in RNA splicing and processing, the same processes enriched in the intersection sets (Figure 5A and Supplementary Figure 11). Inspecting the heatmap and Revigo results (as well as Supplementary Figure 11 and Supplementary Data 5–7) show that, in addition to the RNA splicing and processing terms, other processes were frequently enriched in the sub-modules in the same green cluster, such as proteolysis, protein catabolism, autophagy, and ribosomal biogenesis.

Additionally, the grey cluster displayed a considerable overlap in sub-modules with the green cluster and contained terms related to DNA repair and DNA damage response. This suggests that these processes change in a coordinated fashion during aging, which was already indicated with the enrichment of the 427-set (Figure 3).

These terms’ sub-modules are grouped in the blue, lilac, and yellow columns/submodule clusters (Figure 5C).

Subsequently, we counted the number of DEGs and DCGs in the intersection sets in each module in each tissue (Figure 5D). Most genes in the intersection sets are grouped in the first (and largest) modules, which contain sub-modules in Figure 5C.

Only the heart and skin display different behaviors. This is expected in the case of heart since the sub-modules were defined using only DEGs and DCGs with FDR < 0.1, which were very few.

Only M7 submodules related to skin are found in the blue cluster (Figure 5C). Almost all M2 submodules, which concentrate more genes in the intersection sets, had no enriched terms (Supplementary Data 5–7), explaining why they were not clustered in Figure 5A.

M1 and M5 skin submodules are present in Figure 5C but are contained in the yellow submodule cluster, which, although it has a pattern of terms similar to the blue submodule cluster, has a higher number of terms from the gray cluster, related to DNA repair and response to DNA damage.

The M3 skin module has enriched terms in the submodules for positive DEGs, negative DCGs, and negative DECGs (Supplementary Data 5–7). However, they are more enriched in terms of ribosomal biogenesis, translation regulation, and RNA processing. Only in the submodule of negative DECGs do terms related to RNA splicing appear enriched at FDR < 0.1.

These results show that although the genes in the intersection sets were the only ones changing across tissues, they associate with modules that contain other genes that are changing in a tissue-specific manner, related to similar processes, which act in a coordinated fashion with the genes in the intersection sets.

Protein-protein interaction prediction network analysis

Co-expression relationships do not identify the molecular mechanisms through which biological interactions manifest. Therefore, to determine whether there are physical associations between the protein products of these genes, we constructed protein-protein interaction (PPI) networks using STRING. The STRING database has curated and predicted protein-protein interactions based on several inferences.

Several proteins in the analysis have strong evidence of being part of physical complexes (Figure 6A, 6B). Many of them cluster into very well-defined physical complexes. As expected from our previous results, several interact in the spliceosome.

Additionally, many ribosomal proteins (both cytosolic and mitochondrial) are present, along with translation initiation and elongation factors. Several proteins involved in transcription and chromatin remodeling, as well as proteins involved in ubiquitination and proteolysis, are clustered. The PP1 complex highlighted is also involved in regulating RNA splicing [33].

The complexes identified in this analysis delineate a continuous path, with the transcription of mRNAs being followed by RNA splicing and processing and subsequent translation in the ribosomes, after which proteins are sent to the ER. During this process, damaged or misfolded mRNAs and proteins are degraded. The degradation of proteins is represented in the network by proteins that promote the ubiquitination of other proteins and tag them for degradation, as well as some protease-related proteins. In mitochondria, in addition to the mitochondrial ribosome, some proteins are components of the electron transport chain. Interestingly, only components of Complexes I (NADH dehydrogenase) and III (cytochrome c oxidase) seem to be altered. Similar results are seen in the whole network (Supplementary Figure 12).