Comprehensive transcriptional landscape of porcine cardiac and skeletal muscles reveals differences of aging

Aging significantly affects the cardiac muscle (CM) and skeletal muscles (SM). Since the aging process of CM and SM may be different, high throughput RNA sequencing was performed using CM and SM in different age conditions to evaluate the expression profiles of messenger RNA (mRNA), long non-coding RNA (lncRNA), micro RNA (miRNA), and circular (circRNA). Several mRNAs, lncRNAs, and miRNAs were highly expressed and consistently appeared in both ages in one of the two muscle tissues. Gene ontology (GO) annotation described that these genes were required for maintaining normal biological functions of CM and SM tissues. Furthermore, 26 mRNAs, 4 lncRNAs, 22 miRNAs, and 26 circRNAs were differentially expressed during cardiac muscle aging. Moreover, 81 mRNAs, 5 lncRNAs, 79 miRNAs, and 62 circRNAs were differentially expressed during aging of skeletal muscle. When comparing the expression profiles of CM and SM during aging, the senescence process in CM and SM was found to be fundamentally different. In addition, we assessed multi-group cooperative control relationships and constructed circRNA-miRNA-mRNA co-expression networks in muscular aging. In conclusion, our findings will contribute to the understanding of muscular aging and provide a foundation for future studies on the molecular mechanisms underlying muscular aging.

associated with neural functions, whereas up-regulation of genes was associated with inflammatory responses and heat shock factors [2]. Similarly, the epigenetics of aging may alter the mechanisms involved in gene expression, causing inherited DNA sequences and noncoding RNAs to contribute in defining the broad variety of aging phenotypes [7]. Several non-coding RNAs, including microRNAs (miRNAs) and long non-coding RNAs (lncRNAs), have recently gained attention because of their vital functions [8,9]. Likewise, in some studies, the regulatory mechanisms of miRNAs [10][11][12] and lncRNAs [13,14] were reported in aging mammals. Moreover, many studies documented the expression of circular RNAs (circRNAs) in both humans and mice, which were previously thought to be barely expressed [15][16][17]. At the cellular level, it has been confirmed that circRNAs are involved in chronic CD28-associated CD8(+)T cell aging [18]. Moreover, circRNAs were systematically expressed in skeletal muscle and were differentially expressed in advanced aged muscle, suggesting that circRNAs may affect muscle function [19]. In conclusion, circRNAs may play an important role in the regulation of aging.
Presently, one in every nine individuals worldwide is aged 60 years or older and this number is projected to increase to one in five by the year 2050. This reflects that the number of elderly people will be higher than that in the past [20], and highlights that the challenges of age-related diseases are on the rise. For example, heart failure has been reported as a serious threat to the elderly population and will continue to increase in both primary and secondary care [21]. Furthermore, slowing of processing speed is a general hallmark of human aging, which is related to the decrease in mass, strength, and contraction rate of skeletal muscle during aging [22]. The sarcopenia is primarily responsible for this decline, which becomes increasingly prevalent with age [23]. The cardiac muscle (CM) and skeletal muscle (SM) are the most important muscle tissues that perform vital body functions or support to the human body. The aging process negatively affects the ability and efficiency of CM and SM, and a decrease in mass, strength, rate of contraction and sarcopenia in SM, cardio myopathy, and cardiac fibrosis are all associated with aging [22]. Therefore, these diseases may be closely related with CM and SM. CM and SM are derived from the mesoderm layer of embryonic germ cells, which differentiate into different types of cells and have different biological functions [24], which may be result of different mechanisms of molecular regulation. Indeed, when comparing the transcriptional profiles of CM cells with that of SM cells, differences in gene expression were observed [25]. This demonstrated that aging of CM and SM may be fundamentally different, and suggested that several genes including that of mRNAs and non-coding RNAs may be involved in aging of CM and SM, whereas other mRNAs and non-coding RNAs may play a role in different age-related diseases in these muscle tissues. To our knowledge, only few studies have focused on changes of the transcriptome during CM or SM aging [26,27], which showed that core mechanisms of cardiac aging were conserved from flies to mammals. Moreover, changes in mitochondrial metabolism and protein handling, and contractile function were conserved in hearts of Drosophila, thereby mirroring several aging processes in the mammalian hearts [26], and indicate that SM aging-related miRNAs may contribute to muscle aging with a positive regulation of transcription, and metabolic processes and kinase activity in mouse [27]. Although, these studies have expanded our knowledge on the role of age-related genes in the aging process, comprehensive changes of the transcriptome during aging of CM or SM were not studied, and no attention was paid to the differences between aging in CM and SM. Therefore, the expression profiles of mRNAs, lncRNAs, miRNAs, and circRNAs during aging of CM and SM need to be better understood, the mechanism of mutual regulation and the age-related diseases with these RNAs need to be explored, and the different regulatory mechanisms of the comprehensive transcriptional landscape in the process of aging between CM and SM needs to be revealed.
In previous studies, it was reported that pigs have a genomic structure that, except for anatomic, physiologic and biochemical similarities [28], is very similar to that of humans and mice [29], and has been found to be a valuable model for studying human diseases [30,31]. However, only a few aging studies have been performed using pigs as an animal model. Approximately 1 years of a pig's life is equal to 5 years of human life, resulting in an average life expectancy of the pig of 15-20 years. Therefore, pigs can serve as a relevant model for studying the process of aging because of their longer lifespan compared to that of rodents and similar metabolic features, cardiovascular system, and proportional organ sizes relative to that of humans [32][33][34].
To better understand the transcriptional regulatory mechanism of the aging muscle, we used pigs as our animal model and compared CM and SM in older pigs at the aging stage to the CM and SM in younger pigs at maturation stage. High-throughput sequencing was used to perform an integrated multi transcriptome-wide profiling (mRNA, miRNA, lncRNA, and circRNA) analysis. This procedure enabled us to identify several differentially expressed genes or regulatory elements and unclear multi-group cooperative control networks in the aging muscle. Taken together, these data may help explain muscle age-associated changes in the transcriptional pattern during aging.

rEsuLts changes in the phenotypic traits of cardiac and skeletal muscle tissues during aging
To evaluate changes in phenotypic traits of CM and SM tissues during aging, we collected CM and SM Oncotarget 1526 www.impactjournals.com/oncotarget from two groups of sows (8 years old at aging stage and 180 day old young sows at maturation stage, n = 2 per group), that were prepared for hematoxylin-eosin staining (Supplementary Figure 1). Based on a longitudinal section of CM (10 μm), we observed that the space between muscle fibers of aged individuals was wider compared to that of younger individual (Supplementary Figure  1A). Next, we compared the area of SM fibers of young individual and aged individual, and found that the area of SM fibers of an aged individual (66.50 μm 2 ) was smaller compared to that of a young individual (68.39 μm 2 ) (Supplementary Figure 1B).

transcriptome profile of cardiac and skeletal muscles
To assess to what extent aging induces changes at the transcriptional level, we collected CM and SM from the two groups and performed transcriptomewide profiling, including mRNA, miRNA, lncRNA, and circRNA via high-throughput sequencing. For RNA sequencing libraries, an average of ~72.14 million clean reads were obtained from each of the eight samples and 65.5%~71.8 % of these reads were uniquely aligned to the reference genome Ensemble Susscrofa 10.2 (Supplementary Table 1). Moreover, for small RNA sequencing libraries, approximately 9.25~17.00 million clean reads were obtained from each of the eight samples and 72.7%~90.3% of these reads were uniquely aligned to the reference genome Ensemble Susscrofa 10.2 (Supplementary Table 2).
In the samples, a total of 20339 mRNAs were detected ( Figure 1E). In addition, we identified 7090 and 4402 circRNAs in CM and SM, respectively (Supplementary Dataset 1).
Based on the expression of the transcripts, eight samples, including old cardiac muscle (OCM) and old skeletal muscle (OSM), young cardiac muscle (YCM) and young skeletal muscle (YSM) were separated ( Figure  1A), indicating tissue-specific expression patterns at the transcriptional levels. The Principal Component Analysis (PCA) showed that the mRNAs and lncRNAs transcripts from OCM or OSM had a high degree of variance between each biological replicate, whereas the YCM and YSM results were clustered ( Figure 1B). PCA analyses for miRNA expression resulted in similar findings (Supplementary Figure 2A). Furthermore, we used unsupervised Euclidean matrix plots for expressed mRNA and lncRNA transcripts to determine the variance between muscle of old and young individuals ( Figure 1C). The data showed that a similar tissue (CM or SM) had a closer distance than same age. Moreover, a closer distance value was found between young individuals compared to that of older individuals. A similar result was observed based on unsupervised Euclidean matrix plots of miRNAs (Supplementary Figure 2B). Thus, the results indicated that the transcriptional levels of OCM and OSM were more variable than that of their younger counterparts.

clustering special profiles of gene expression
To further determine the pathways which involved in the development of the CM and SM transcriptome, we used the k-means clustering algorithm to systematically screen highly expressed RNAs, which were consistently present in both ages and identified the Gene Ontology (GO) terms with the help of the GO consortium database (described in the Methods section). After grouping 105 mRNAs in the CM-cluster 1 (Figure 2  , which basic function was oxygen transport as one of the globin chain components of hemoglobin A [35]. Mutant HBB caused sickle cell anemia [36], and had a higher expression level in OCM compared to YCM, suggesting that the aging CM needs to increase oxygen transport to maintain normal levels of oxidation. Furthermore, we identified the PNPLA2 (ENSSSCG00000012841) gene that was related to positive regulation of triglyceride catabolic processes, and was significantly up-regulated (log 2 FC = 3.03390467391754) during aging CM. Studies have shown that mutations in PNPLA2 caused neutral lipid storage disease combined with myopathy (NLSDM) [37][38][39][40][41]. Down-regulated genes were related to the structural molecule, gene transcript and calcium binding capacity, such as structural molecule activity (GO: 0005198), aging (GO: 0007568), RNA polymerase II activity (GO: 0001055), and calcium ion binding (GO: 0005509) ( Figure 3A, Supplementary Dataset 6). Downregulation of the gene CDH9 (ENSSSCG00000020845, log 2 FC = -9.034974243) (Supplementary Dataset 4), which is a classic cadherin [42], was enriched for calcium ion binding. It encodes neuronal cell-adhesion [43] and had a lower expression level in OCM compared to YCM, indicating it may be related to the autonomic nervous system (controlling cardiac function), which decreases with age [44]. In addition, we observed hyaline degeneration in OCM by HE-stained muscle sections ( Figure 3E). This reflected the structural changes of CM tissue during aging, which may be related to inflammaging [45].
Similarly, we obtained 81 DEGs in OSM vs YSM (OYSM), which showed 51 up-regulated genes and 30 down-regulated genes ( Figure 3B, Supplementary Dataset 4). GO analysis demonstrated that these upregulated genes were enriched for immune response, cell differentiation, metabolic process, glucose transport, and apoptotic process, such as the defense response to virus (GO: 0051607), skeletal muscle cell differentiation (GO: 0035914), lipid catabolic process (GO: 0016042), oxygen metabolic processes (GO: 0072592), glucose transport (GO: 0015758), and negative regulation of apoptotic process (GO: 0043066) ( Figure 3D, Supplementary Dataset 6). We found that the genes enriched for "defense response to virus", such as ISG15 (log 2 FC = 6.71813789075041), IFIT3 (log 2 FC = 5.37430318712509), IFIT1 (log 2 FC = 5.05205229449539) (Supplementary Dataset 4), showed higher expression levels in OSM compared to YSM. Previous studies have shown that these genes were closely related to innate immunity responses [46,47]. Moreover, over activation of ISG15 was found to uniquely associated with autoimmune disease dermatomyositis [48]. In addition, we observed inflammatory cell infiltration in HE-stained sections of OSM ( Figure 3F), which may confirm that the changes in gene expression were related to "defense response to virus" tissue. Down-regulated genes were related to structural proteins, DNA binding, neuron development, response to stress and smoothened binding, such as caveolae (GO: 0005901), sequence-specific DNA binding (GO: 0043565), cerebellar neuron development (GO: 0098749), regulation of response to osmotic stress (GO: 0047484) and smoothened binding (GO: 0005119) ( Figure 3D, Supplementary Dataset 6). For example, the gene HOXA9 (ENSSSCG00000028997) (log 2 FC = -9.96363818460901) (Supplementary Dataset 4) enriched for sequence-specific DNA binding, was significantly down-regulated in SM, which confirmed that the deletion of HOXA9 improved satellite cell function and muscle regeneration in aged mice [49].  Figure 4D). While the target genes of down-regulated SM-DELs were closely related to DNA biding and gene transcription, such as sequence-specific DNA binding (GO: 0043565), regulation of transcription, DNA-template (GO: 0006355), some target genes were involved in cerebellar neuron development (GO: 0098749), and relaxation of muscle (GO: 0090075) ( Figure 4E, Supplementary Dataset 7).
As the initial exploration of functional implications of differences in lncRNA across genomes, we investigated whether lncRNAs co-regulated the differences in target gene expression during muscle aging. Based on the results of DELs and DEGs, we investigated 7 co-expression DELs and their differentially expressed co-target genes in both tissues (Table 1). In CM, the HBB (DEG) was a co-target of the two up-regulated DELs (MSTRG.51934 and MSTRG.86430). HBB was related to oxygen transport and caused sickle cell anemia by a mutation in HBB [35,36]. In SM, the two up-regulated DELs (MSTRG.212099 and MSTRG.49389) and their target genes were coupregulated. ARRDC2 is a target of MSTRG.212099 and GO analysis showed that ARRDC2 was involved in "cytoplasmic vesicle". A previous study showed that ARRDC2 was dynamically distributed throughout the plasma membrane and endocytic system, and may play a role in cargo protein trafficking within the endocytic system by mediating the discrete G protein coupled receptor (GPCR) [50]. Moreover, IFIT3 GO enriched for "defense response to virus", which was related to the innate immunity response, showed a maximum of Pearson correlation in four target genes of MSTRG.49389. In SM, the 3 down-regulated lncRNAs (MSTRG.82735, MSTRG.59589 and MSTRG.8641) demonstrated four differentially expressed co-target genes, including NOS1, OTP, NANOS1, and ENSSSCG00000026009. Based to GO enrichment analysis, these co-target genes were related to oxidoreductase activity, neurohypophysis development, cerebellar neuron development, and receptor-mediated endocytosis. Several studies showed that NOS1 gene polymorphism was associated with asthma [51], OTP played an essential role in the specification of neuronal cell lineages in the developing hypothalamus [52] and NANOS1 was required for maintaining oocyte production [53].   Oncotarget 1531 www.impactjournals.com/oncotarget the formation of cardiac fibrosis, and that expression of miR-29b-3p significantly decreased with an increase in fibrosis marker and collagen protein. Moreover, in cardiac fibroblasts, Na/K-ATPase signaling regulates collagen synthesis through miR-29b-3p [54]. Here, we showed that miR-29b-3p was up-regulated during aging of CM and may contribute to cardiac fibrosis in the aging heart [22,55]. Interestingly, we demonstrated that both miR-493-5p and its target FGB ( Figure 5B) were significantly and differentially expressed in SM, showing that miR-493-5p (log 2 FC = -4.977874815) was down-regulated while FGB (ENSSSCG00000008997, log 2 FC = 8.954205569) was  Dataset 4). FGB is the beta component of fibrinogen, and when vascular injury occurs, fibrinogen is cleaved by thrombin to form fibrin, which is the most abundant component of blood clots [56].
We also identified 26 DECs in CM and of these DECs, 20 circRNAs were up-regulated while 6 circRNAs were down-regulated (Supplementary Figure 3C Oncotarget 1533 www.impactjournals.com/oncotarget and circRNA005658) were simultaneously expressed in both CM and SM (Supplementary Dataset 4). The circRNA020853 was up-regulated in CM but downregulated in SM, while circRNA005658 was downregulated in CM but up-regulated in SM. In addition, circRNA020877 was shown to be down-regulated in both CM and SM.

construction of the circrnA-mirnA-mrnA coexpression network
It has been shown that circRNAs may act as competing endogenous RNAs (ceRNAs), which regulate the function of miRNAs [57,58], suggesting that circRNAs and their target miRNAs may be co-expressed in the aging muscle. The co-expression network between circRNAs and their target miRNAs were predicted by the Pearson correlation coefficient and based on DECs and DEMs data [59]. For CM, a network map was constructed, which contained 9 DEMs, 4 DECs, and 11 relationships ( Figure 6A, Supplementary Dataset 8). Within the network, all miRNAs were co-related with simultaneous up-regulation of circRNAs; circRNA005698, with 7 miRNAs, showing most relations. Similarly, a complex network map was constructed for SM, which contained 26 DEMs, 17 DECs, and 69 relationships ( Figure 6B, Supplementary Dataset 8). In the network, several upor down-regulated miRNAs were co-related with up-or down-regulated circRNAs, and circRNA014844, with 11 miRNAs, showing the most relations.
These findings prompted us to further analyze the circRNA-miRNA-mRNA co-expression network. In the network, we found that FGB could be targeted by miR-493-5p, suggesting that it may be a crucial factor mediated by circRNAs-miR-493-5p axes. According to the DECs, DEMs, and DEGs data, we created the circRNA (circRNA014844, circRNA011308, and circRNA018401)-miR-493-5p-FGB co-expression pathway in SM (Figure 7). In the network, the down-regulated miR-493-5p was positively co-related with three down-regulated circRNAs (circRNA014844, circRNA011308, and circRNA018401), while it was negatively co-related with the up-regulation of its target gene FGB.

the expression of genes validated by q-Pcr
Here, we performed q-PCR to validate the high throughput sequencing results. We selected 9 differentially Oncotarget 1534 www.impactjournals.com/oncotarget expressed genes to performed q-PCR, including 3 mRNAs, 2 lncRNAs, 2 miRNAs, and 2 circRNAs (Supplementary Figure 4). Expression of several selected genes was significantly different during aging, which may be closely related to biological functions or function in a network, such as PNPLA2, HOXA9, FGB (Supplementary Figure 4A), lncRNA MSTRG.8641 (Supplementary Figure  4B), miRNA miR-493-5p (Supplementary Figure 4C), and circRNA circRNA014844 (Supplementary Figure 4D). Similarly, several genes were randomly selected from the differentially expressed genes, including MSTRG.143042 (Supplementary Figure 4B), miR-29b-3p (Supplementary Figure 4C), and circRNA010340 (Supplementary Figure  4D). The data indicated that the expression patterns were highly consistent between the two methods ( Supplementary Figure 4), for example, according to the high throughput sequencing, miR-29b-3p was up-regulated in both aging CM and SM, and the expression patterns were highly consistent with that of the q-PCR results. The q-PCR results confirmed the high reproducibility and reliability of the gene expression profiling in our study.

dIscussIon
To our knowledge, this is the first time that the profiles of transcriptome expression were comprehensively analyzed and the differences compared between CM and SM during aging. In our study, we found significant differences in gene expression for CM and SM. Several genes were highly expressed in one tissue and hardly expressed in the other ( Figure 1A), indicating the different functions at a general transcriptome level in these muscle types. In addition, we identified genes for tissue-specific expression, and found that 4.9% of 20339 genes (n = 994) were expressed in CM and 4.4% (n = 886) in SM ( Figure 1D). Furthermore, we analyzed these genes and identified that several highly expressed genes did not show an obvious change during aging. Indeed, GO annotation showed these genes were required for maintaining normal biological function of CM and SM. This provides new information through which we can view the differences in functions of CM and SM at the mRNA, lncRNA, and miRNA levels.
In our study, we demonstrated an increased individual variation by analyzing mRNAs and lncRNAs ( Figure 1B) in CM and SM during aging. Similar findings were found in aging mouse liver by analyzing mRNAs [60]. Additionally, we examined miRNAs to verify this variation (Supplementary Figure 2A) and similar results were obtained. Moreover, this was the first time to report variations of lncRNA and small RNA, indicating that increased individual variation during aging was a common feature that not only reflected multiple organizations but also reflected multi-layer of the transcriptome.
We identified several age-related genes for CM and SM during aging and found more up-regulated genes compared to down-regulated genes. These results were in line with similar findings reported in previous studies [26,[60][61][62]. Perhaps, when the body ages the ability of several metabolic functions decreased and could be compensated by increasing associated genes expression. During aging, we found a higher number of differentially expressed genes in SM compared to CM, which may be related to fact that SM performs more complex biological functions than CM. We found that differentially expressed genes were closely related to the biological function of tissue and certain diseases. Here, we found that the PNPLA2 (ENSSSCG00000012841) gene was related to positive regulation of triglyceride catabolic process and was significantly up-regulated (log 2 FC = 3.03390467391754) during CM aging. Previous studies reported that mutations in PNPLA2 could cause neutral lipid storage disease with myopathy (NLSDM) [37][38][39][40][41]. This implied that PNPLA2 expression increased during aging of the heart, and improved triglyceride metabolism which may be important for preventing NLSDM in the process of biological aging. We found that HOXA9, which is related to sequencespecific DNA binding ( Figure 3D, Supplementary Dataset 6), was significantly down-regulation ((log 2 FC = -9.96363818460901)) in SM during aging, and was negligibly expressed in older muscle tissue. In a previous study, it was demonstrated that the inhibition of aberrant chromatin activation or deletion of HOXA9 improved the satellite cell, which can differentiate and merge to form new mature multi-nucleated muscle cells [63], function and muscle regeneration in aged mice [49]. Likewise, HOXA9 significantly declined in SM, which may indicate that the SM needs to increase the rate of muscle regeneration to repair damaged muscle fibers during aging. However, the capability of muscle reconstruction decreases at a later age [64].
In our study, we found that miR-493-5p and the target of miR-493-5p decreased and increased, respectively during SM aging. Further, according to the co-expression network of circRNA-miRNA, we found that circRNA014844, circRNA011308, and circRNA018401 were closely related to miR-493-5p expression. Therefore, we constructed an mRNA-miRNA-circRNA interaction network, which showed that FGB is related to procoagulation [56]. Therefore, we speculated that as the muscle tissue ages, procoagulation in the muscle may become weak and more fibrinogen needs to be produced to accelerate the generation of blood clots. FGB, miR-493-5p, and circRNAs may play a key role in regulating procoagulation during aging. Similarly, studies have shown that mutations in the FGB gene could lead to afibrinogenemia [65][66][67], and FGB tyrosine nitration is a prothrombotic risk factor [68]. A higher thrombotic risk in elderly was partially due to a higher procoagulation state [69]. Thus, the increased FGB during aging found in our study may be contributed to this risk.

Animals
Four healthy female Yana pigs (an indigenous Chinese pig) were used in this study. These pigs included two 180 day-old young sows and two 8 year-old sows. There were no direct and collateral blood relationships within the last 3 generations among the 4 pigs. The piglets were weaned at 28 ± 1 day of age. A starter diet was administered from day 30 to day 60 after weaning, and included 3. 40 Mcal kg −1 metabolizable energy with 20.0% crude protein (11.5 g/kg lysine). From day 61 to day 120, the pigs were fed a diet, containing 14.0 MJ/kg of metabolizable energy consisting of 18.0% crude protein (9.0 g/kg lysine). From day 121, pigs received a diet containing 13.5 MJ/kg of metabolizable energy and 16.0% crude protein (8.0 g/kg lysine). Animals were allowed access to food and water ad libitum and were maintained under the same conditions. The night before slaughtering, pigs were fastened and given two hours for rest after transportation. Then, pigs were stunned electrically at 90 V and 50 Hz for 10 s, and exsanguinated as essential to ameliorate pain.

samples preparation
All samples used in this study were collected according to the guidelines for care and use of experimental animals established by the Ministry of Science and Technology of China. From each pig, the following tissues were collected: CM, longissimus dorsi muscle (LDM), SM. Tissues were rapidly separated from each carcass, immediately frozen in liquid nitrogen, and stored at -80°C until RNA extraction. Total RNA was extracted using TRIzol Reagent (Invitrogen, CA, USA), and treated with DNase and purified using an RNeasy Mini Kit (Qiagen, Valenica, CA, U.S.A.). The concentration and quality of RNA was assessed on the Agilent Bioanalyzer 2100 system.

rnA sequencing and data analysis
A total of 5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the rRNA-depleted RNA by NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina ® ((New England Biolabs Ipswich, Massachusetts, USA) following the manufacturer's guidelines and library quality was assessed using an Agilent Bioanalyzer 2100 system. After cluster generation, libraries were sequenced on an Illumina HiSeq 4000, and 150 bp paired-end reads were generated. Clean reads were obtained after removing reads containing the adapter, reads containing ploy-N and low quality reads from raw data. Clean reads were aligned to the Ensemble (Susscrofa 10.2) using TopHat2 (v2.0.14) [70] with default parameters.
StringTie software [71] was used to assemble the mapped reads of each sample, which at least existed in one of both replicates. Transcripts included BLASTed (evalue = 1e-10) to Ensembl, and mapped transcripts were directly described as known lncRNA or mRNA. Here, we used Salmon (v0.6.0) [72] to calculate transcripts per million (TPMs) of both lncRNAs and coding genes in each sample. Then, Coding Potential Calculator (CPC) (0.9) [73] and Pfam Scan (v1.5) [74] were performed to analyze the coding potential of transcripts. Transcripts predicted to have coding potential by either/all of the tow tools were filtered out, and those without coding potential served as our candidate set of novel lncRNAs.

small rnA sequencing and data analysis
A total of 5 μg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® ((New England Biolabs Ipswich, Massachusetts, USA) following the manufacturer's guidelines, and index codes were added to attribute sequences to each sample. Library quality was assessed using an Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina, Inc, California) according to the manufacturer's guidelines. After cluster generation, the library preparations were sequenced on an Illumina MiSeq platform, and 50 bp single-end reads were generated.
The miRBase21 was used as reference, and software mirdeep2 [75] was used to obtain the potential miRNA, and to predict novel miRNA. In detail, we used miRanda (v3.3a) software [76] with default parameters and cutoffs (Score S ≥ 140 and Energy E ≤-20.0) to predict miRNA target. Subsequently, miRNA expression levels were estimated by TPMs by the following criteria: normalization formula: normalized expression = mapped read count/total reads*1000000.

clustering profiles of gene expression
In this study, the k-means clustering algorithm was applied to analyze the expressed genes. For each cluster [77], gene normalized expression level was presented in log 2 (tpm+1).

differential expression analysis
Differentially expressed mRNAs, lncRNAs, miRNAs and circRNAs were discovered using the edgeR package [78] running in the R programming environment. The edgeR user guide was followed as detailed in the "classic analysis" section. A differential expression q-value < 0.05 and fold change > 2 assigned as differentially expressed in two different comparisons www.impactjournals.com/oncotarget (OCM vs. YCM, OSM vs. YSM). For DECs, we used a standard p-value < 0.05 and fold change > 2, fold change was log2 transformed ratio of average expression between two groups, calculated by edgeR.

Gene ontology enrichment analysis
GO (Gene Ontology) enrichment analysis was implemented by the GOseq R package, GO terms with corrected p-value < 0.05 were considered significantly enriched by differentially expressed genes.

circrnA-mirnA co-expression network
The circRNA-miRNA co-expression network was constructed based on the correlation analysis between the differentially expressed circRNA and miRNAs. The expressions of differentially expressed circRNAs and miRNAs were analyzed by Pearson's correlation coefficient. The absolute coefficient value of 0.8 between a circRNA and a miRNA was considered relevant for network construction. A p-value of < 0.05 was considered statistically significant.

circrnA-mirnA-mrnA co-expression network
The mRNA-miRNA interactions were predicted by miRanda (v3.3a), Minimum Free Energy (MFE) threshold values for miRanda was set to -20 kcal/mol. We first analyzed the results of mRNA-miRNA and circRNA-miRNA co-expression network, then constructed the circRNA-miRNA-mRNA co-expression network.

Quantitative Pcr validation
cDNA was synthesized using the oligo (dT) and random 6-mer primers provided in the PrimeScript RT Master Mix kit (TaKaRa, Shiga, Japan). q-PCR was performed using the SYBR Premix Ex Taq kit (TaKaRa, Shiga, Japan) on a CFX96 Real-Time PCR detection system (Bio-Rad Laboratories, Inc., Hercules, CA, U.S.A.). The q-PCR validation was carried out using three biological replicates. The primer pairs used are presented in Supplementary Table 3. Three endogenous control genes (porcine GAPDH, ACTB, and U6 snRNA) were used in this assay. The ΔΔCt method was used to determine the expression level of objective mRNAs, miRNAs, lncRNAs, and circ RNAs.

concLusIons
In summary, in the present study the known expression profiles of transcriptome for CM and SM are explored during aging. We grouped several tissuespecific high expression genes at both ages, and analysis results indicated that transcriptional levels in active gene regulation manifests itself and altered transcriptome profiles in a tissue-specific manner. Furthermore, we found that several differentially expressed genes were closely related to some serious health conditions and that may several genes may play a crucial role in regulating the functions during aging. We also compared the expression profiles of CM and SM during aging and found significant differences between the two muscle types, indicating the senescence process of CM and SM is fundamentally different. In addition, we constructed the circRNA-miRNA-mRNA co-expression networks in muscular aging cycle. Our findings contribute to the understanding of muscular aging and provide a foundation for future studies of the molecular mechanisms underlying muscular aging.