Muscle fiber type-dependence effect of exercise on genomic networks in aged mice models

Principal component analysis (PCA) plots and heatmaps

To understand the expression profiles of the identified genes, PCA plots and heatmaps are presented. SOL and GAS were divided into two clusters, as illustrated by the dendrogram (Figure 1A). DEGs were mostly obtained in the comparison between types of muscle fibers (SOL vs. GAS) (Figure 1A and Supplementary Tables 1, 2). PCA plots also indicated that two clusters (SOL vs. GAS) were clearly separated by principal component 1 (PC1, 69.2%) (Figure 1B). Considerable separation between the young and old groups was achieved by principal component 2 (PC2, 5.5%) (Figure 1C).

Identification of DEGs between YC vs. OC in SOL and GAS and commonly expressed DEGs

To investigate the effect of aging on SOL and GAS, DEG analysis was performed. DEGs were determined following >1.5-fold change in random expectation (p < 0.05). A total of 286 DEGs were obtained with aging in both SOL and GAS, including 196 DEGs in SOL OC/YC and 119 DEGs in GAS OC/YC (Supplementary Table 1). Of the 196 DEGs in SOL OC/YC, 90 genes were upregulated, and 106 were downregulated (Supplementary Table 1). However, among 119 DEGs in GAS OC/YC, 71 genes were upregulated, and 48 were downregulated (Supplementary Table 1). For each of SOL and GAS, the 30 DEGs with the largest absolute value of fold change are plotted (Figure 2A, 2B).

To determine how many of the genes were changed owing to aging in both SOL and GAS, we identified 29 commonly expressed DEGs (Figure 2C and Table 1). Of the 29 common DEGs, only cluster of differentiation 74 (Cd74) and TYRO protein tyrosine kinase-binding protein (Tyrobp) exhibited a different pattern between SOL and GAS (Table 1). For Cd74 and Tyrobp, the fold change of OC/YC showed downregulation in SOL, whereas it showed upregulation in GAS (downregulation in SOL OC/YC and upregulation in GAS OC/YC) (Table 1). Excluding Cd74 and Tyrobp, 27 of 29 overlapped DEGs showed similar patterns between SOL and GAS (Table 1). For seven common DEGs, the fold change of OC/YC exhibited upregulation in both SOL and GAS (upregulated in both SOL OC/YC and GAS OC/YC) (Table 1). In contrast, for 20 of 29 common DEGs, the fold change of OC/YC displayed downregulation in both SOL and GAS (downregulation in both SOL OC/YC and GAS OC/YC) (Table 1).

Identification of DEGs between YC vs. YE in SOL and GAS and commonly expressed DEGs

To elucidate the effect of exercise on SOL and GAS in the young groups, DEG analysis was performed. DEGs were identified following >1.5-fold change in random expectation (p < 0.05). A total of 42 DEGs were obtained as a result of exercise in young groups, including 31 DEGs in SOL YE/YC and 13 DEGs in GAS YE/YC (Table 2). Of the 31 DEGs in SOL YE/YC, 11 genes were upregulated, and 20 were downregulated, whereas among 13 DEGs in GAS YE/YC, seven genes were upregulated, and six were downregulated (Table 2). For each of SOL and GAS, the 30 DEGs with largest absolute value of fold change are plotted (Figure 3A, 3B and Table 2).

To identify how many of the genes were commonly changed owing to exercise in young groups between SOL and GAS, we checked for the overlapping DEGs in a Venn diagram. The number of DEGs commonly expressed was two (Figure 3C). Of the two common DEGs, nuclear receptor subfamily 1 group D member 1 (Nr1d1) exhibited a similar pattern between SOL and GAS, meaning the fold change of YE/YC showed downregulation in both SOL and GAS (downregulation in both SOL YE/YC and GAS YE/YC) (Table 2). In contrast, phosphoinositide-3-kinase regulatory subunit 1 (Pik3r1) showed a different pattern between SOL and GAS. Namely, for Pik3r1, the fold change of YE/YC exhibited upregulation in SOL, whereas it exhibited downregulation in GAS (upregulation in SOL YE/YC and downregulation in GAS YE/YC) (Table 2).

Identification of DEGs between OC vs. OE in SOL and GAS and commonly expressed DEGs

To elucidate the effect of exercise on SOL and GAS in the old groups, DEG analysis was performed. DEGs were identified following a >1.5-fold change in random expectation (p < 0.05). A total of 30 DEGs were obtained by exercise in the old groups, including 20 DEGs in SOL OE/OC and 13 DEGs in GAS OE/OC (Table 3). Out of 20 DEGs in SOL OE/OC, 10 genes were upregulated, and 10 were downregulated, whereas among 13 DEGs in GAS OE/OC, six genes were upregulated, and seven were downregulated (Table 3). For each of SOL and GAS, the 30 DEGs with largest absolute value of fold change are plotted (Figure 4A, 4B).

To determine how many of the genes were commonly altered by exercise in old groups between SOL and GAS, we explored commonly expressed DEGs between SOL and GAS (DEGs of SOL OE/OC and GAS OE/OC). Three common DEGs were found (Figure 4C), of which patatin-like phospholipase domain-containing 3 (Pnpla3) exhibited a similar pattern between SOL and GAS (Table 3). Specifically, for Pnpla3 the fold change of OE/OC displayed upregulation in SOL and GAS (upregulation in both SOL OE/OC and GAS OE/OC) (Table 3). However, the hemoglobin subunit beta-1 (Hbb-b1) and chromatin target of PRMT1 (Chtop) showed different patterns (Table 3). For Hbb-b1, the fold change of OE/OC exhibited upregulation in SOL, whereas it exhibited downregulation in GAS (upregulation in SOL OE/OC and downregulation in GAS OE/OC) (Table 3). In the case of Chtop, the fold change of OE/OC displayed downregulation in SOL but upregulation in GAS (downregulation in SOL OE/OC and upregulation in GAS OE/OC) (Table 3).

Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using DEGs of SOL/GAS in young and old groups

To clarify the differences between SOL and GAS by utilizing DEGs in SOL YC/GAS YC and SOL OC/GAS OC, GO term and KEGG pathway enrichment analyses were performed. The number of DEGs in SOL YC/GAS YC was 1060, while in SOL OC/GAS OC, it was 1052; these DEGs were used for enrichment analyses (Supplementary Table 2). For the top 10 GO terms by SOL YC/GAS YC, the pathways associated with fatty acid metabolism were significantly enriched, including fatty acid beta-oxidation (GO:0006635), fatty acid metabolic process (GO:0006631), and lipid metabolic process (GO:0006629) (Figure 5A and Table 4). For the top 10 KEGG pathways in the young groups, the enrichment scores of fatty acid degradation (mmu00071) and fatty acid metabolism (mmu01212) were significant (Figure 5B and Table 5). Regarding valine, leucine, and isoleucine degradation (mmu00280), this pathway had only 14 upregulated DEGs in SOL YC/GAS YC (Table 5). In addition, out of 19 DEGs of glycolysis/gluconeogenesis (mmu00010), 15 genes were downregulated, and four were upregulated (Table 5).

For the top 10 GO terms by SOL OC/GAS OC, the gene set related to lipid metabolism, such as fatty acid metabolic process (GO:0006631) and lipid metabolic process (GO:0006629), was significantly enriched (Figure 5C and Table 6). Regarding the top 10 KEGG pathways by SOL OC/GAS OC, the enrichment scores of fatty acid metabolism (mmu01212) and fatty acid degradation (mmu00071) were significant (Figure 5D and Table 7). For pathways associated with glucose metabolism, three gene sets belonged to the top 10 KEGG pathways in the old groups (Figure 5D and Table 7). These pathways included glucagon signaling (mmu04922), glycolysis/gluconeogenesis (mmu00010), and insulin signaling (mmu04910) (Figure 5D and Table 7).

Gene set enrichment analysis (GSEA) for apoptosis in SOL and GAS

To investigate the association between aging and exercise and apoptosis in SOL and GAS, we analyzed data at the gene set levels. For the comparison between young and old groups, HALLMARK_APOPTOSIS pathway was more enriched in GAS OC than in GAS YC (Figure 6A and Table 8), whereas in SOL, the enrichment score (ES) between SOL YC and SOL OC was not significant (ES = 0.36, p = 0.099). For the comparison between sedentary control and exercise in old groups, the HALLMARK_APOPTOSIS pathway was more enriched in GAS OC than in GAS OE (Figure 6B and Table 9). In SOL, the enrichment score between SOL OC and SOL OE was not significant (ES = −0.25, p = 0.615). In the young groups, there were no significant difference between the control and exercise groups in both SOL and GAS.

GSEA analysis for myogenesis in SOL and GAS

To elucidate the relationship between aging and exercise and myogenesis in SOL and GAS, we conducted GSEA. For the comparison between young and old groups, HALLMARK_MYOGENESIS pathway was more enriched in SOL YC than in SOL OC (Figure 7A and Table 10), whereas in GAS, the enrichment score between GAS YC and GAS OC was not significant (ES = −0.25, p = 0.955). For the results of GSEA using sedentary control and exercise in young groups, the HALLMARK_MYOGENESIS pathway was more enriched in GAS YC than in GAS YE (Figure 7B and Table 11). Similar to GAS, this pathway was also more enriched in SOL YC than in SOL YE (Figure 7C and Table 12). For the comparison between control and exercise in old groups, HALLMARK_MYOGENESIS pathway was more enriched in GAS OE than in GAS OC (Figure 7D and Table 9), whereas in SOL, it was more enriched in OC groups than in OE groups (Figure 7E and Table 12).

Effect of muscle fiber type on gene expression is greater than that of aging and exercise

PCA plots and heatmaps indicate that the separation by muscle fiber type is most conspicuous among the clusters by muscle fiber type, aging, and exercise. For separation by fiber type, SOL was divided from GAS by PC1 (69.2% of total variance) (Figure 1B). Prior studies performing PCA confirmed that the gene expression of fast-twitch muscles is different from that of slow-twitch muscles by PC1 (46% of total variance) [27], which our findings are compatible with. However, although both EDL and psoas are mainly composed of fast-twitch fibers, EDL exhibited clear separation from psoas by PC2, accounting for 27% of the total variance [27]. This result suggests that in addition to muscle fiber type distribution, other factors may have a significant effect on gene expression. Since the whole muscle lysates used in this study may include blood vessels, connective tissues, and nerves in addition to muscle fibers, EDL may display different gene expression patterns from psoas owing to differences in these factors. Therefore, to elucidate the effects of muscle fiber type on gene expression more clearly, microgenomic profiles, such as single-cell sample analysis, may have to be performed [28]. After the muscle fiber type, clusters by aging were the most prominent. The young and old groups were noticeably separated by PC2 (5.5%) (Figure 1C). In addition, previous data have shown that metabolomics datasets are remarkably clustered into young and old groups [29], which our results agree with.

In the current study, the exercise groups were not clearly divided (Figure 1B, 1C). According to previous research performing PCA, microarray data were prominently separated into exercise and sedentary groups [30]. Presumably, these conflicting results between our study and the one performed previously may stem from the differences in the period of exercise intervention and tissues used in gene expression analysis. For the period of exercise, our protocol required 4 weeks, while exercise intervention in the prior study lasted for 7 weeks. To determine whether exercise for 4 weeks is long enough to elicit the differences in gene expression profiles, further research should be conducted. Furthermore, the experiment in our study relied on skeletal muscles such as SOL and GAS, whereas the left ventricle was utilized in the previous research [30]. To identify whether the cardiac muscle is more responsive than SOL and GAS when intervening with the same exercise regimen, a follow-up study is needed.

Common DEGs with aging are related to apoptosis, circadian rhythm, and immune signaling

We performed DEG analysis in SOL and GAS to elucidate the differences due to aging. There were 29 common DEGs that were associated with aging (Figure 2C), of which cell death-inducing DNA fragmentation factor alpha like effector A (Cidea) was upregulated in both SOL and GAS with age (Table 1). Cidea has been confirmed to be involved in apoptosis and lipid metabolism [31–33]. In the skeletal muscle of rats, apoptotic signaling was enhanced, and Cidea expression was significantly increased after burn injury [31]. For lipid metabolism, Cidea-knockout mice have a lean phenotype, whereas mice with induced transgenic expression of Cidea exhibit a healthy obese phenotype [32, 33]. Although it has been confirmed that Cidea expression is decreased by aging in adipose tissue [34], it is not known whether Cidea expression changes with aging in skeletal muscle. To clarify what influence the alteration of Cidea expression in skeletal muscles has on phenotype, further studies are necessary. In addition, collagen genes (Col1a1, Col1a2, Col3a1, Col4a2, Col5a1, Col5a2, Col5a3, and Col6a1) were downregulated in both SOL and GAS with aging. Prior studies also observed that the concentration of collagen was decreased by aging [35, 36].

D-box-binding PAR BZIP transcription factor (Dbp) and Nr1d1 (also known as Rev-erb-alpha) were downregulated in both SOL and GAS and are referred to as “clock-controlled genes” (Table 1). Considering that reduced Dbp expression is correlated with the early aging phenotype in mice [37], there seems to be an interaction between Dbp and aging. Furthermore, overexpression of Nr1d1 increases mitochondrial function and exercise capacity [38]. In contrast, the deficiency of Nr1d1 leads to the upregulation of atrophy-related genes and increases skeletal muscle atrophy [39]. As we have confirmed here that Nr1d1 expression is decreased in skeletal muscles owing to aging (Table 1), further investigation is needed to scrutinize the interaction between Nr1d1 and muscle atrophy induced by aging. However, circadian genes exhibit differences in expression depending on the timing of measurement, tissues, and species [40]; therefore, it is necessary to consider these factors as well.

Cd74 and Tyrobp (also known as Dap12) display different expression patterns between SOL and GAS by the effect of aging (Table 1). Importantly, these genes are involved in immune functions. Migration inhibitory factor (the primary ligand of Cd74) knockout mice reported to have lived markedly longer than control mice [41]. Tyrobp is mainly expressed on innate immune cells, including monocytes, neutrophils, dendritic cells (DCs), and natural killer (NK) cells [42]. Because Cd74 and Tyrobp were barely addressed in terms of aging and exercise effects, it is meaningful to determine whether the differences in their expression between SOL and GAS elicit dissimilar characteristics of immune function.

Common DEGs with exercise are associated with circadian rhythm and insulin signaling in young groups and hemoglobin subunit in old groups

In young groups, there were two DEGs commonly expressed owing to exercise: Nr1d1 and Pi3kr1 (Table 2). In our study, Nr1d1 expression was downregulated in both SOL and GAS (Table 2). However, as prior research has shown that acute exercise increases Nr1d1 expression [43], our findings seem to contradict these previous data. Putatively, these inconsistent results may be due to differences in the exercise period (acute vs. chronic). Considering that Nr1d1 has a positive correlation with mitochondrial function and oxidative capacity in skeletal muscles, the increase in Nr1d1 expression by acute exercise may be due to the need to temporarily improve mitochondrial function. On the other hand, the reason chronic exercise increases Nr1d1 expression may be because chronic exercise increases oxidative capacity so the exercise group can possibly maintain homeostasis without the influence of Nr1d1. Nonetheless, the correlation between Nr1d1 and exercise period remains to be explored. Pi3kr1 displayed different expression patterns in SOL and GAS as a result of exercise (Table 2). For instance, Pi3kr1 expression was upregulated in SOL and downregulated in GAS (Table 2). According to a previous study, Pi3kr1 is indispensable for insulin signaling, and the selective deletion of Pi3kr1 augmented insulin sensitivity in fast-twitch muscle (EDL) but not in SOL [44]. As the effect of exercise on glucose uptake differs depending on muscle fiber type composition [45], further research is warranted to determine whether the differences in Pi3kr1 expression partially contribute to the differences in insulin resistance according to muscle fiber type.

In old groups with exercise, common DEGs included Pnpla3 and Hbb-b1 (Table 3). Pnpla3 expression was upregulated in both SOL and GAS (Table 3). However, as Pnpla3 is mainly located in the liver and retina [46], the alteration of Pnpla3 expression in skeletal muscle does not seem to be important. Hbb-b1, on the other hand, showed different patterns between SOL and GAS (Table 3). With exercise, Hbb-b1 expression was upregulated in SOL and downregulated in GAS (Table 3). Hbb-b1 is linked to the formation of hemoglobin, and the mutation of Hbb-b1 is associated with abnormal oxygen transportation [47]. As Hbb-b1 is involved in oxygen transportation and exhibits the most prominent differences between SOL and GAS, further studies should be conducted to clarify the interaction between this protein and muscle fiber type.

Differences in metabolism according to muscle fiber type and significant upregulation of BCAA degradation in SOL compared to that in GAS

To identify the differences between SOL and GAS, we performed GO analysis using DEGs between these two muscle types. For GO and KEGG pathways in all age groups, fatty acid metabolism was upregulated in SOL, while gene sets associated with glucose metabolism were downregulated (Figure 5A–5D and Tables 4–7). Considering that the enzyme activities of fatty acid metabolism are higher in slow-twitch muscle than in fast-twitch muscle and that glucose metabolism is invigorated in fast-twitch muscle [48], our findings are compatible with those shown previously. For the KEGG pathway in young groups, the valine, leucine, and isoleucine degradation pathways were significantly upregulated in SOL (Figure 5B and Table 5). Considering that branched-chain amino acid (BCAA) catabolism leads to a decrease in myofiber size [49], there is a need to study whether the difference in BCAA degradation elicits different patterns of atrophy based on muscle fiber type.

Necessity for different interventions depending on the muscle fiber type due to different patterns of apoptosis and myogenesis with aging

In fast-twitch muscle, our results showed that the apoptotic pathway was significantly enriched in the old group compared to the young group (Figure 6A and Table 8), whereas no significant difference was observed with aging in terms of myogenesis. Considering that the muscle mass-to-body weight ratio is decreased in fast-twitch muscle with age [50], the cause of atrophy may be largely owing to changes in the apoptotic process due to aging, not myogenesis. Unlike GAS, no significant difference in the apoptotic process due to aging was seen in SOL, while myogenesis was significantly enriched in the young group compared to the old group (Figure 7A). According to a prior study, the cross-sectional area of slow-twitch muscles, such as SOL, is decreased by aging in rats [51]. Given the results of our study, atrophy of slow-twitch muscle may primarily result from the alteration of myogenesis with age and not from the change in apoptosis. As the cause of atrophy by aging seems to differ depending on muscle fiber type composition, appropriate intervention according to fiber type should be considered. As the interaction between aging-induced atrophy and muscle fiber type has rarely been investigated regarding both apoptosis and myogenesis, further research is needed to clarify the differences in these two processes by aging depending on fiber type at the cellular level.

Mitigation of apoptosis and augmentation of myogenesis in old groups by chronic aerobic exercise in GAS but not in SOL

We found that exercise intervention for 4 weeks had different effects on the apoptosis and myogenesis depending on the muscle fiber type. In old groups, the apoptotic pathway was significantly enriched in the sedentary control compared to the exercise group in GAS (Figure 6B and Table 9), whereas no significant difference was seen in SOL. Our results are consistent with the results of previous research, which show that exercise significantly attenuates apoptosis in fast-twitch muscle but not in slow-twitch [52]. According to prior study, age-related muscle atrophy mainly relies on intrinsic apoptotic pathway and this process are different between slow and fast-twitch muscle [51, 53–56]. In detail, caspase-9, which is a key molecule in the intrinsic apoptotic pathway, were significantly up-regulated in fast-twitch muscle by aging, whereas were not in slow-twitch muscle [50]. Since it has been not determined whether chronic exercise results in down-regulation of apoptosis by modulating the amount and activity of caspase-9 in fast-twitch muscle, further study need to be performed.

In the young groups, there was no significant difference due to exercise in either SOL or GAS. Previous research has shown that chronic exercise induces different apoptotic processes depending on the muscle fiber type in 4-week-old rats [57]. Putatively, the difference in our results with those found previously may be due to the differences in the type of exercise (concentric vs. eccentric). The exercise protocol in the prior study relied on eccentric exercise by decreasing the slope of the treadmill [57], while our protocol calls for a 6° incline. As eccentric exercise elicits more muscle damage than concentric exercise [58], previous studies using eccentric exercise may have exhibited more conspicuous differences in the apoptotic process depending on muscle fiber type than our current study. Whether these conflicting results stem from the differences in the type of exercise remains to be investigated.

In the old groups, myogenesis was more enriched in exercise groups than in the sedentary control in GAS (Figure 7D and Table 9), whereas there were no significant differences seen in SOL. Our results are consistent with those shown previously where 8 weeks of treadmill exercise significantly increase the level of the myogenin, transcription factor involved in myogenesis, in fast-twitch muscle but not in slow-twitch muscle in aged rats [59]. Considering that myogenesis is involved in muscle regeneration [60], a follow-up study should be performed to identify whether chronic treadmill exercise enhances the regenerative capacity of fast-twitch muscles more so than slow-twitch muscles in old groups.

In the young groups, myogenesis was not more enriched in the exercise group than in the sedentary group in both SOL and GAS (Figure 7B, 7C, Tables 11 and 13). According to a previous study, treadmill exercise for 13 weeks doesn’t induce to significant increase in myogenic clones of GAS in young group, whereas myogenic clones was significantly elevated in old group [61]. Therefore, aerobic exercise can mitigate the decrease in myogenesis by aging, but does not seem to up-regulate myogenic pathway more actively in young group. Since there are not many studies done in this regard, further studies are needed to clarify. Taken together, future related studies should be mindful that the interaction of muscle fiber type with aging and exercise is relevant with respect to apoptosis and myogenesis.

Animals

All procedures were based on previously published data [63]. Briefly, the 9-week-old C57BL/6 young male mice and 84-week-old mice were purchased from Central Lab. The groups were divided into young control (YC), young exercise (YE), old control (OC), and old exercise (OE) groups. Mice were individually housed in standard conditions with food and water ad libitum. Each group comprised three mice and all mice were subjected to RNA sequencing analysis. Details on animal care were provided in a previous study [63].

Exercise protocol

Before 4 weeks of treadmill exercise, adaptation exercise was conducted once a day for 3 days. Mice were familiarized with the treadmill for 15 min/session at a 0 m/min for 3 min, 5 m/min for 2 min and 8 m/min for 10 min and 6° incline. In 4 weeks of treadmill exercise, mice allocated to perform treadmill running were subject to 6° incline and for warm up, 2 min speed of 0 m/min, speed of 5 m/min, 8 m/min, 10 m/min for 1 min each, and then 12 m/min for 30 min at first week, 2 m/min were increased every week, and cool down at 5 m/min for 2 min for 1 session (37 min). 2 session/day were performed. Between the session, there were at least 1 h break time was given.

SOL and GAS RNA extraction

All procedures were based on a previous study [63]. Briefly, total RNA was extracted from SOL and GAS of mice using a total RNA extraction kit (Ribozol, Amresco) following the manufacturer’s protocol. The remaining processes, including cDNA synthesis, check for purity and integrity were presented in previously published data [63].

RNA-Seq library preparation and sequencing

Total RNA concentration was calculated by Quant-IT RiboGreen (Invitrogen, #R11490). To assess the integrity of the total RNA, samples are run on the TapeStation RNA screentape (Agilent, #5067–5576). Only high-quality RNA preparations, with RIN greater than 7.0, were used for RNA library construction. A library was independently prepared with 1ug of total RNA for each sample by Illumina TruSeq Stranded mRNA Sample Prep Kit (Illumina, Inc., San Diego, CA, USA, #RS-122-2101). The first step in the workflow involves purifying the poly-A containing mRNA molecules using poly-T-attached magnetic beads. Following purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments are copied into first strand cDNA using SuperScript II reverse transcriptase (Invitrogen, #18064014) and random primers. This is followed by second strand cDNA synthesis using DNA Polymerase I, RNase H and dUTP. These cDNA fragments then go through an end repair process, the addition of a single ‘A’ base, and then ligation of the adapters. The products are then purified and enriched with PCR to create the final cDNA library. The libraries were quantified using KAPA Library Quantification kits for Illumina sequencing platforms according to the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, #KK4854) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies, #5067–5582). Indexed libraries were then submitted to an Illumina NovaSeq (Illumina, Inc., San Diego, CA, USA), and the paired-end (2 × 100 bp) sequencing was performed by the Macrogen Incorporated.

Data pre-processing and quality check

Raw data (FPKM) were processed using the Bioconductor preprocessCore package with default parameters for filtering, signal +1 and logarithm transformation, and quantile normalization [64]. DEG analysis was performed using fold change calculation and the Bioconductor genefilter package [65].

Bioinformatics tools