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

Skeletal muscles are made up of various muscle fiber type including slow and fast-twitch fibers. Because each muscle fiber has its own physiological characteristics, the effects of aging and exercise vary depending on the type of muscle fiber. We used bioinformatics screening techniques such as differentially expressed gene analysis, gene ontology analysis and gene set enrichment analysis, to try to understand the genetic differences between muscle fiber types. The experiment and gene expression profiling in this study used the soleus (SOL, slow-twitch muscle) and gastrocnemius (GAS, fast-twitch muscle). According to our findings, fatty acid metabolism is significantly up-regulated in SOL compared to GAS, whereas the glucose metabolism pathway is significantly down-regulated in SOL compared to GAS. Furthermore, apoptosis and myogenesis patterns differ between SOL and GAS. SOL did not show differences in apoptosis due to the aging effect, but apoptosis in GAS was significantly up-regulated with age. Apoptosis in GAS of old groups is significantly reduced after 4 weeks of aerobic exercise, but no such finding was found in SOL. In terms of myogenesis, exercise intervention up-regulated this process in GAS of old groups but not in SOL. Taken together, muscle fiber type significantly interacts with aging and exercise. Despite the importance of the interaction between these factors, large-scale gene expression data has rarely been studied. We hope to contribute to a better understanding of the relationship between muscle fiber type, aging and exercise at the molecular level.


INTRODUCTION
Type I (slow-twitch) and type II (fast-twitch) muscle fibers make up the mammalian skeletal muscle, which is a very heterogeneous tissue [1]. Different factors, such as aging, exercise, and diseases, affect type I and type II fiber differently. Previous studies have shown that aging has a different impact depending on the muscle fiber type composition. Type II fibers significantly decreased with age, whereas type I fibers were unaffected [2][3][4]. Furthermore, glycolysis and glycogen metabolism were also up-regulated in type I fibers as they aged, whereas they were down-regulated in type II fibers [2]. As a result, the interaction between aging and muscle fiber type is important in terms of muscle atrophy and metabolism, though the exact cellular mechanism is still unknown.
One of the most effective strategies for rejuvenation, cognitive function, and disease prevention, such as metabolic syndrome, is exercise [5][6][7]. Exercise has different effects depending on muscle fiber type, in addition to aging. Acute exercise significantly increased the expression of lipolytic complex genes (ATGL, AGING G0S2) in the soleus (SOL) and red section of the gastrocnemius (GAS), whereas it decreased the expression of lipolytic complex genes (ATGL, CGI- 58, G0S2) in the white section of the gastrocnemius (GAS) [8]. The effect of exercise on muscle hypertrophy varies depending on the muscle fiber type. Resistance training improves the muscle mass-to-body ratio in fast-twitch muscles but not in slow-twitch muscles [9]. Moreover, exercise causes changes in the composition of muscle fiber types. Type IIA fibers were found to be increased in previous studies, while type I fibers remained unchanged or decreased [10][11][12]. Endurance training, on the other hand, can cause a shift toward type I fibers [13,14]. To summarize, exercise has a close relationship with muscle fiber type, as evidenced by numerous studies, and the relationship between exercise and muscle fiber type should be investigated to better understand the numerous benefits of exercise.
Among skeletal muscles, SOL and GAS have a different muscle fiber type composition than the rest of the skeletal muscles. In mouse SOL contains a high percentage of slow twitch fibers (type I 37.42 ± 8.20%, IIA 38.62 ± 6.81%, IIAD 18.74 ± 6.95% and IID 5.69 ± 3.09%), whereas GAS is primarily made up of fast twitch fibers (type IIB 54.42 ± 8.11%, IIDB 19.37 ± 2.98%, IID 2.26 ± 2.24%, IIAD 12.40 ± 2.34%, IIA 5.73 ± 3.24% and I 5.74 ± 2.55%) [15]. Because SOL and GAS have remarkably different muscle fiber type composition, SOL and GAS exhibit different patterns with aging. In rats, SOL was not significantly changed in muscle weights by aging, while GAS was decreased [16]. In humans, SOL thickness did not differ in older vs. younger groups, but GAS atrophied in older vs. younger groups [17]. SOL has the characteristics of a slow-twitch muscle, whereas GAS has the characteristics of a fast-twitch muscle, due to the fact that fast-twitch fiber tends to atrophy more with age than slow-twitch fiber.
In addition to the effect of aging, exercise has different effects on SOL and GAS. SOL has more DEGs induced by exercise than GAS in terms of number. The stressresponse gene activating transcription factor 3 (Atf3) [18], was only up-regulated in SOL and not in GAS [19]. Although strenuous exercise causes an increase in apoptosis in both SOL and GAS, the patterns of apoptosis are different (intrinsic vs. extrinsic pathway) [20]. Furthermore, the effects of exercise on brainderived neurotrophic factor (BDNF, known to promote neuron survival) and sequestosome 1 (SQSTM1, a receptor that accumulates when autophagy is inhibited) differ between SOL and GAS [21]. Short-term exercise increased BDNF expression in SOL but not in GAS (for 5 days). SQSTM1 gene expression was increased in SOL, but not in GAS, despite chronic exercise (for 8 weeks) [22]. In summary, exercise has a different effect on between SOL and GAS at the molecular level.
Taken together, muscle fiber type, aging and exercise relevantly interact with each other. Because the effects of aging and exercise differ significantly depending on muscle fiber type, and muscle fiber type composition is linked to diseases such as obesity, diabetes, and cancer [23][24][25][26], it's critical to look into the interaction between muscle fiber type and various factors. Despite these requirements, large-scale gene expression data has rarely been examined in terms of the effects of aging and exercise on different muscle fiber types. As a result, we used bioinformatics tools to perform systematic gene expression profiling in an attempt to elucidate the genetic relationship between muscle fiber type, aging, and exercise.

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).  Genes with * are different pattern between SOL and GAS. Genes with * are different pattern between SOL and GAS. AGING 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).  Genes with * are different pattern between SOL and GAS.

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.  Table 3.

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).

DISCUSSION
Muscle fiber type interacts with aging and exercise; aging has a major impact on type II fibers in terms of atrophy, and exercise induces hypertrophy to a different extent depending on fiber type [2,9]. As such, many researchers have scrutinized the interaction of muscle fiber type with aging and exercise. However, large sets of RNA-sequencing data have rarely been analyzed to clarify the relationship between these factors. Therefore, by utilizing such data and bioinformatics, we attempted in the current study to elucidate how the muscle fiber type interacts with aging and exercise at the genetic level.

Enriched in SOL OC
There was no geneset that satisfied the following criteria (p value < 0.05, FDR < 0.25) Size, the number of corresponding genes to gene sets. Abbreviations: ES: Enrichment Score; NES: Normalized Enrichment Score.

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][32][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 atrophyrelated 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][54][55][56]. In detail, caspase-9, which is a key molecule in the intrinsic apoptotic pathway, were significantly upregulated 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 downregulation 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 slowtwitch 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.

CONCLUSION
There is a myriad of conflicting results from studies on exercise and aging. As research scientists, we are obligated to identify the causes of inconsistent outcomes in scientific experiments. In this study, we demonstrated that gene expression patterns differ significantly according to muscle fiber type (slowtwitch vs. fast-twitch), making it imperative that this be considered when investigating exercise and aging. Furthermore, to contribute to taking a step closer to a unanimous conclusion, elucidating how exercise frequency, intensity, time, type, and intervention period interact with muscle fiber type distribution would be invaluable.
This study has limitation. SOL and GAS exhibits different activation patterns with increasing demands of force and speed [62]. In detail, as exercise intensity increased, the EMG amplitude ratios of SOL to GAS decreased [62]. Thus, in order to exclude mechanical differences, exercise is performed at an intensity that induced the same level of activation of SOL and GAS.
The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE198266 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GS E198266).

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].

PCA analysis, heatmap, and Venn diagram
PCA was performed using the built-in R function prcomp() [66]. The Github ggbiplot and ggplot2 packages were used to plot a PCA graph [67]. Preprocessed data were used for PCA. The Bioconductor ComplexHeatmap package was utilized to create a heatmap [68]. For the heatmap in Figure 1, preprocessed data were screened based on the following criteria: After the median value and standard deviation (SD) for each gene were calculated using the log2 expression value, genes with a median value of <5 and SD of the bottom 25% were excluded from the analysis. Of the 18,795 total genes, 1268 passed this screening. For the heatmap representing the results of GSEA, the data were quantile normalized without performing signal + 1 and logarithm transformation. The Venny 2.1 online service was used to create a Venn diagram [69].

DEG analysis
DEGs were identified using the Bioconductor gene filter package in R [65]. To test whether there was a statistically significant difference between the groups, an independent t-test was applied to each gene of the pre-processed data. The statistical threshold for significance was a p-value < 0.05 and fold change >1.5 for DEGs derived from SOL OC/YC, SOL OE/OC, SOL YE/YC, GAS OC/YC, GAS OE/OC, and GAS YE/YC. The DEGs originating from SOL OC/GAS OC, SOL OE/GAS OE, SOL YC/GAS YC, and SOL YE/GAS YE were applied to p-value < 0.05 and fold change >2.0.

GO and KEGG pathway enrichment analysis
Functional annotation analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) and KEGG. The criteria for statistically significant GO terms and KEGG pathways were FDR-corrected p-values < 0.05. To visualize the results of the enrichment analysis, the GOplot in R was utilized [70]. To facilitate interpretation, the z-score was calculated using the following formula:

GSEA
GSEA is a computational method that determines whether a priori defined set of genes shows statistically significant, concordant differences between two biological states. A gene set that satisfied p-value < 0.05 and FDR <0.25 was determined to be significant.