Research Paper Volume 12, Issue 6 pp 5071—5090

Comprehensive analysis of circular RNA profiles in skeletal muscles of aging mice and after aerobic exercise intervention

Mingwei Guo 1, 2, *, , Jin Qiu 2, *, , Fei Shen 3, *, , Sainan Wang 2, , Jian Yu 2, , Hui Zuo 2, , Jing Yao 2, , Sainan Xu 2, , Tianhui Hu 2, , Dongmei Wang 2, , Yu Zhao 3, , Yepeng Hu 1, , Feixia Shen 1, , Xinran Ma 1, 2, , Jian Lu 3, , Xuejiang Gu 1, , Lingyan Xu 1, 2, ,

  • 1 Department of Endocrine and Metabolic Diseases, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, Zhejiang, China
  • 2 Shanghai Key Laboratory of Regulatory Biology, Institute of Biomedical Sciences and School of Life Sciences, East China Normal University, Shanghai, China
  • 3 Key Laboratory of Adolescent Health Assessment and Exercise Intervention, Ministry of Education, East China Normal University, Shanghai, China
* Equal contribution

received: July 13, 2019 ; accepted: March 9, 2020 ; published: March 17, 2020 ;

https://doi.org/10.18632/aging.102932
How to Cite

Copyright © 2020 Guo et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Aging induces gradual accumulation of damages in cells and tissues, which leads to physiological dysfunctions. Aging-associated muscle dysfunction is commonly seen in aged population and severely affects their physical activity and life quality, against which aerobic training has been shown to exert antagonizing or alleviating effects. Circular RNAs (circRNAs) play important roles in various physiological processes, yet their involvement in aging-associated muscle dysfunction is not well understood. In this study, we performed comprehensive analysis of circRNAs profiles in quadriceps muscles in sedentary young and aging mice, as well as aging mice with aerobic exercise using RNA sequencing. Our results identified circRNAs altered by factors of aging and aerobic exercise. Their host genes were then predicted and analyzed by gene ontology enrichment analysis. Importantly, we found that circBBS9 featured decreased levels in aging compared to young mice and elevated expression in exercise versus sedentary aging mice. Besides, we performed GO and KEGG analysis on circBBS9 target genes, as well as established the circBBS9-miRNA-mRNAs interaction network. Our results indicate that circBBS9 may play active roles in muscle aging by mediating the benefits of aerobic training intervention, thus may serve as a novel therapeutic target combating aging-associated muscle dysfunction.

Introduction

An increase in life expectancy in modern world brings an ‘aged society’, in which a substantial aging population poses challenges both medically and financially [1, 2]. Thus, more and more research interests are piqued toward the combat against aging. Generally, aging is defined as the age-dependent physiological decline that affects all living organisms. Aging undermines multiple major organs and plays a profound role in the onset of neurodegenerative diseases, cardiovascular diseases, metabolic disorders, as well as a loss in muscle and bone mass [37].

Recently, metabolic fitness emerges as a novel player in the arena of combating aging. For instance, the prevalence of calorie-enriched diets and sedentary lifestyle of today’s society has been shown to be closely linked with the onset and deteriorating of various age-associated diseases [8, 9]. Evidences from genetic modified mice models revealed that alterations in genes controlling metabolic homeostasis have significant impacts on overall longevity of the animals [10, 11]. Last but not least, metabolic interventions, such as caloric restriction, intermittent fasting and exercise, confer multiple health benefits ranging from stronger skeleton muscle and cardiac function to improved metabolic fitness and cognitive functions, and most importantly, increased lifespan in multiple species [12, 13].

Skeletal muscle comprises about 40% of total body mass in mammals and consumes a large proportion of fuel molecules in both resting and active states. These characteristics attribute a central role to skeletal muscle in maintaining whole body physical fitness and metabolic homeostasis [1418]. During aging, skeletal muscle undergoes aging-associated muscle dysfunction, which contributes greatly to the disrupted metabolic homeostasis. Multiple factors affect the process of aging-associated muscle dysfunction. Nowadays, the rapid development and wide application of RNA sequencing technology have enabled an extensive understanding of the role of non-coding RNAs in muscle dysfunctions. For instance, a number of micro RNAs are reported to show significant influence on the physiology and pathology of muscle fibers [1921]. Circular RNA (circRNA), a subfamily of the non-coding RNAs, is ubiquitously expressed in eukaryotes with tissue- and developmental-stage-specific characteristics, which feature a covalently closed loop in its structure. Compared to miRNA, circRNAs are exceptionally stable due to their circular structure that lacking free 5’or 3’ ends, thus endows them superior potential as signaling molecules and biomarkers [22]. Indeed, circRNAs exhibit critical functions under both physiological and pathological scenarios through diverse regulations on gene transcription, pre-mRNA splicing, mRNA translation and protein functions, as well as their classic roles as miRNA sponges [23]. Overall, circRNAs are attractive candidates as biomarkers and therapeutic targets of aging-associated muscle dysfunction [24].

The aim of the study was to find potential circRNA involved in muscle aging and exercise via profiling the expression signature and heterogeneity pattern of circRNAs in quadriceps femoris muscles of sedentary young and aging mice, as well as aging mice with aerobic training. Our results identified circBBS9 and its target gene pathways affected by the advance of age, as well as reversed by aerobic exercise, which may offer novel insights into the biomarkers and the pathogenesis of aging-associated muscle dysfunction.

Results

Description of circRNA profiles in quadriceps femoris muscle from young, aging and aging with aerobic exercise group of mice

To access the extents of muscle dysfunction, six muscle samples per group from mice of young, aging and aging plus aerobic exercise were examined for gene markers indicative of mitochondrial homeostasis and muscle atrophy. Compared to young control mice, we found significant mitigation in mitochondrial program (Pgc1a, Mfn1 and Atpase) and elevation of atrophy markers (Foxo3 and Atrogin) in aging mice, while aging mice with aerobic training showed marked amelioration in these parameters (Figure 1A) as consistent with previous reports [25]. Next, using RNA-seq technique, we evaluated the circRNA profiles of quadriceps femoris muscles in three groups of samples following the workflow shown in Figure 1B. The sequencing identified 4336 circRNAs in total in three groups (Supplementary Table 1). The numbers of circRNA distributed in genome were the highest in chromosome 2 and reached nadir in chromosome Y, with most chromosomes had a synchronized distribution of about 200 counts (Figure 2A). The majority of the identified circRNAs were less than 1500 nucleotide (nt) in length (Figure 2B). The relationship between circRNAs and their coding genes were summarized and classified into three categories: 85.56% were exonic, 5.60% were intronic and 8.84% were others types (Figure 2C). Besides, all circRNAs were evenly located on DNA plus and minus strand (Figure 2D).

Validation of aging and aerobic training model in mice and workflow of circRNA analysis scheme. (A) The Qu muscle expression level of genes in atrophy and mitochondrial functionality among groups of young, aging and aging with aerobic exercise (n=6 per group). (B) Workflow of circRNA analysis scheme.

Figure 1. Validation of aging and aerobic training model in mice and workflow of circRNA analysis scheme. (A) The Qu muscle expression level of genes in atrophy and mitochondrial functionality among groups of young, aging and aging with aerobic exercise (n=6 per group). (B) Workflow of circRNA analysis scheme.

Differences and characterizations of circRNA expression profile. (A) Chromosomal distributions of annotated circRNAs. (B) Predicted spliced length of circRNAs. (C) The circRNA were classified into three types according to the relationship of the genomic loci with their associated coding genes. (D) Distribution of circRNA in sense (+) and antisense (-) strand of DNA.

Figure 2. Differences and characterizations of circRNA expression profile. (A) Chromosomal distributions of annotated circRNAs. (B) Predicted spliced length of circRNAs. (C) The circRNA were classified into three types according to the relationship of the genomic loci with their associated coding genes. (D) Distribution of circRNA in sense (+) and antisense (-) strand of DNA.

Identification of differential circRNAs in quadriceps muscles

We firstly investigated the differential circRNAs in quadriceps muscle between young and aging group, which may be used as biomarkers for aging-associated muscle dysfunction. The volcano plot filtered and identified the differentially expressed circRNAs with statistical significance (Figure 3A). The threshold of exhibiting fold change is 2.0 and p-values below 0.05. Between the young and aging group, 49 circRNAs showed significant differential expression, with 28 circRNAs up-regulated and 21 circRNAs down-regulated (Supplementary Table 2). Meanwhile, comparison between aging mice and aging mice with aerobic treadmill training identified 21 significantly changed circRNAs, among them 10 were up-regulated and 11 were down-regulated (Figure 3B, Supplementary Table 3). The top changed circRNAs were shown in Figure 3C, 3D.

Overview of altered circRNA in different sequencing group. (A) Volcano plots showing differential expression of circRNAs of Aging group compared with Young group. Differentially expressed circRNAs with fold change > 2 and p B) Volcano plots showing differential expression of circRNAs of Aging plus Exercising compared to Aging group. Differentially expressed circRNAs with fold change > 2 and p C) The top 10 upregulated and 10 downregulated circRNA based on the log2 fold change of Aging group compared with Young group. (D) The top 10 upregulated and 10 downregulated circRNA based on the log2 fold change of Aging plus Exercising group compared to Aging group.

Figure 3. Overview of altered circRNA in different sequencing group. (A) Volcano plots showing differential expression of circRNAs of Aging group compared with Young group. Differentially expressed circRNAs with fold change > 2 and p < 0.05 were marked in orange and green dots representing up and down regulation separately. (B) Volcano plots showing differential expression of circRNAs of Aging plus Exercising compared to Aging group. Differentially expressed circRNAs with fold change > 2 and p < 0.05 were marked in orange and green dots representing up and down regulation separately. (C) The top 10 upregulated and 10 downregulated circRNA based on the log2 fold change of Aging group compared with Young group. (D) The top 10 upregulated and 10 downregulated circRNA based on the log2 fold change of Aging plus Exercising group compared to Aging group.

To explore the potential functional roles of the significantly enriched circRNAs and their host genes, we performed Gene Ontology (GO) enrichment analysis on the host genes. GO analysis is divided into three parts as biological process (BP), cell component (CC) and molecular function (MF). We found that many biological processes (BP), i.e. muscle contraction and skeletal muscle thin filament assembly, as well as cell component (CC) including myofibril, muscle myosin complex and myosin filament were altered in aging mice compared to young mice. (Supplementary Figure 1). On the other hand, compared to sedentary aging mice, aging mice with aerobic intervention showed altered biological processes (BP) including muscle filament sliding, neuron remodeling, muscle contraction and cell components (CC) including myofibril, myosin complex, cytoskeleton and muscle myosin complex (Supplementary Figure 2). These data suggested a potential cooperation between circRNAs and their host genes in the regulation of muscle functions during aging and aerobic interventions.

Validation of circBBS9 expression in skeletal muscle

Based on the sequencing data, we wonder whether we could define circRNA as potential target molecule that mediates the benefits of aerobic training during aging process by highlighting candidates that featured altered expression in aging while reversed by exercise. Hence, we overlapped circRNAs of three sets of sequencing data, looking for those featured opposite expression patterns between young and aging, as well as between sedentary aging and aging with exercise intervention. Importantly, we found three circRNAs that met our criteria, circBBS9, circATP9b and circMYH8 (Supplementary Table 4).

Next, we verified circBBS9, circATP9b and circMYH8 with qRT-PCR. By qRT-PCR using primers across the branch site [26], we found that muscle circBBS9 level was down-regulated in aging versus young sedentary mice and was restored after aerobic training in aging mice, which was the best match for RNA-seq data (Figure 4A), while circMYH8 and circATP9b levels failed to match the sequencing results (Supplementary Figure 3A). Thus we focused on circBBS9 for detailed analysis.

Verification of the expression of CircBBS9. (A) qRT-PCR verification of the expression of circBBS9 among groups. (B) Schematic diagram of primer design of circBBS9. (C) Identification of circBBS9 in Qu muscle by PCR amplification. (D) Sanger sequencing to verify the amplified products of circBBS9. Data are presented as mean ± SEM and *P

Figure 4. Verification of the expression of CircBBS9. (A) qRT-PCR verification of the expression of circBBS9 among groups. (B) Schematic diagram of primer design of circBBS9. (C) Identification of circBBS9 in Qu muscle by PCR amplification. (D) Sanger sequencing to verify the amplified products of circBBS9. Data are presented as mean ± SEM and *P<0.05, **P<0.01.

According to UCSC genome browser, circBBS9 was spliced from a 416bp pre-mRNA containing two exons and one intron. CircBBS9 has a length of 243bp derived from the joint of two adjacent exons during spicing (Figure 4B). To ensure assay facticity in circRNA, we used RiboNuclease R for full degradation of preponderant linear RNAs and Y-structure RNAs to preserve the circular RNAs. Next, convergent primers designed specifically against the adjacent ends of two exons (interspaced by the intron) and divergent primers against the far end of two exons were tested in RT-PCR analysis in complementary DNA (cDNA) sample and genomic DNA (gDNA) sample from quadriceps muscles. As shown in Figure 4C, convergent primers yielded PCR products of 339bp in gDNA and 166bp in cDNA sample, indicating the exclude of intron during splicing. Meanwhile, unlike in cDNA sample, divergent primers failed to amplify PCR products in gDNA sample, suggesting circBBS9 is a circular RNA derived from the splicing of pre-mRNA. This was confirmed via Sanger sequencing, which showed the existence of a splicing event on AG-GT site (Figure 4D). Interestingly, circBBS9 levels were not altered in brain and heart, two other critical organs underwent senescence and dysfunction during aging (Supplementary Figure 3B). These data suggest that although aging is a systematic process, circBBS9 might be a reliable and specific biomarker for skeletal muscle sarcopenia and aerobic exercise intervention.

Construction of circBBS9-miRNA-mRNA Network and validation of circBBS9 predicted target mRNAs

CircRNAs classically function as miRNA sponges to exert their regulatory effects on gene expression [26]. We thus constructed the circBBS9-miRNA-mRNA network to better unravel its role in muscle dysfunction. The Miranda and RNAhybrid bioinformatics tools were utilized to predict the sponge miRNAs for circBBS9, which predicted 10 miRNAs as its targets (mmu-miR-3100-5p, mmu-miR-6930-3p, mmu-miR-7020-5p, mmu-miR-423-3p, mmu-miR-7079-5p, mmu-miR-383-3p, mmu-miR-6911-5p, mmu-miR-3065-3p, mmu-miR-7028-5p, mmu-miR-7662-5p). Next, using database Targetscan, mRNAs predicted as the targets of at least three miRNAs were overlapped and considered as mRNA targets, which rendered 1558 target mRNAs. These genes were analyzed with GO and KEGG analysis to annotate and speculate their potential functions.

GO analysis on the target genes of circBBS9 highlighted various biological processes, for instance, transcription regulation and protein phosphorylation in BP analysis and protein homodimerization and ATP binding in MF analysis (Figure 5A5C). Subsequent KEGG analysis of target genes emphasized metabolic pathways, PI3K-Akt signaling pathway, MAPK signaling among top 10 enriched pathways (Figure 5D). Of note, metabolic pathways had the highest target gene counts among these signaling pathways, suggesting a critical impact of metabolic changes on muscle functions during aging and exercise intervention. Based on these analyses, the comprehensive atlas of miRNA-mRNA network for circBBS9 was then constructed using Cytoscape software (3.6.1). In detail, we created miRNA-mRNA network if the mRNA were regulated by a greater number (>=3) of miRNAs and had a p-value <0.05. The circBBS9 miRNA-mRNA network was comprised of 7,283 edges between 10 miRNAs and 1,558 mRNAs (Supplementary Table 5, Figure 6A). The degree of distribution of the nodes followed the power-law distribution with a slope of -75.189 and an R-squared value of 0.9798, suggesting the network displayed scale-free characteristics typical of a biological network rather than a random system (Figure 6B). In Figure 6C, compared to mRNAs (4.678), The average degree of the miRNAs was 728.3. Besides, miRNA nodes showed more betweenness centrality compared to mRNA nodes as in Figure 6D. These results suggested that miRNAs have more interactions with other nodes than mRNAs and although miRNAs are small RNAs, they exhibit more topological properties than mRNAs in the network. Interestingly, a few miRNAs acted as major hubs linking other mRNAs, such as mmu-miR-7662-5p, mmu-miR-7028-5p and mmu-miR-6911-5p (degree=940), which were largely beyond the maximum degree node among mRNAs (Hic2, Fbxo41, Dnmt3a, degree=10).

In silico analysis of predicted target genes of circBBS9. (A–C) GO analysis of predicted target genes with top 10 differ gene counts. The horizontal axis is the gene counts for the GO terms, and the vertical axis is the GO terms. (D) KEGG pathway analysis of predicted target gene with top 10 differ gene counts. Selection counts represent the number of entities of target genes directly associated with the listed Pathway.

Figure 5. In silico analysis of predicted target genes of circBBS9. (A–C) GO analysis of predicted target genes with top 10 differ gene counts. The horizontal axis is the gene counts for the GO terms, and the vertical axis is the GO terms. (D) KEGG pathway analysis of predicted target gene with top 10 differ gene counts. Selection counts represent the number of entities of target genes directly associated with the listed Pathway.

CircBBS9 target miRNA-mRNA network analysis. (A) Network of the circRNA-miRNA which have greater interaction score with their target mRNAs. The pink nodes represented miRNA and the green nodes represented mRNAs. (B) Degree distribution of circBBS9 related miRNA-mRNA network. (C, D) The degree and betweenness centrality of mRNAs and miRNAs.

Figure 6. CircBBS9 target miRNA-mRNA network analysis. (A) Network of the circRNA-miRNA which have greater interaction score with their target mRNAs. The pink nodes represented miRNA and the green nodes represented mRNAs. (B) Degree distribution of circBBS9 related miRNA-mRNA network. (C, D) The degree and betweenness centrality of mRNAs and miRNAs.

In order to confirm the network of circBBS9-miRNA-mRNA, we overexpressed circBBS9 in differentiated C2C12 myotubes and examined a few predicted target mRNAs, including Dnmt3a, Dad1, Gys1, Cacnale, Adcy1, Adcyap1r1 and Ctnnd1, which belong to different top predicted pathways. Importantly, we found that circBBS9 overexpression increased these gene expressions (Figure 7A, 7B), in addition to decreased muscle atrophic genes (Foxo3 and Atrogin) and elevated functional mitochondrial genes (Pgc1α, Mfn1 and Atpase) (Figure 7C). Overall, these data suggested that circBBS9 overexpression improve muscle functionality gene programs and may be involved in muscle aging process.

mRNA alternations upon circBBS9 overexpression in differentiated C2C12 myotubes. (A) Expression levels of circBBS9 in differentiated C2C12 myotubes after infection of lentiviral delivery of circBBS9. (B, C) Expression levels of predicted circBBS9 target genes (B) and general muscle functionality genes (C) in differentiated C2C12 myotubes after infection of lentiviral delivery of circBBS9. Data are presented as mean ± SEM and **P

Figure 7. mRNA alternations upon circBBS9 overexpression in differentiated C2C12 myotubes. (A) Expression levels of circBBS9 in differentiated C2C12 myotubes after infection of lentiviral delivery of circBBS9. (B, C) Expression levels of predicted circBBS9 target genes (B) and general muscle functionality genes (C) in differentiated C2C12 myotubes after infection of lentiviral delivery of circBBS9. Data are presented as mean ± SEM and **P<0.01.

Discussion

In this study, we identified and verified circRNAs that showed differential expression between young and old. Analysis on their target genes revealed metabolic pathways among the top affected pathways. We also highlighted circRNAs that featured declined expression in aging compared to young and were reversed by aerobic training. Among them, circBBS9 was verified and its circRNA-miRNA-mRNA regulatory network established. Taken together, our results revealed multiple circRNAs as potential biomarkers of aging and unveiled a previously unappreciated role of circBBS9 in mediating the beneficial effects of aerobic exercise, at lease in aged subjects. Future studies are warranted to study the potential function of circBBS9 as a novel therapeutic target for aging-associated muscle dysfunction.

Aging-associated muscle dysfunction severely affects the health states and life qualities of the elders. As the basic module of skeletal muscle, myocytes, also known as muscle fibers, are packed with mitochondrial to maintain their high ATP demand during exercise for energy sustenance. It is thus vital for myocytes to maintain strict mitochondrial homeostasis for optimal muscle functionality. Notably, muscle mitochondrial homeostasis undergoes dynamic regulation in both. normal and pathological conditions. For example, aged population suffers from aging-associated decrease in muscle mass and strength, with prominent mitochondrial dysfunctions in myocytes. On the other hand, aerobic training is one of the most effective ways in promoting muscle health and mitochondrial functionality by activating key metabolic regulators such as SIRT1, AMPK and PGC1α [17]. However, implementing aerobic exercise in the aged population demands significant caution in frequency and intensity so as to balance its desiring effects with potential damages like bone fractures and cardiovascular events [18]. One alternative in promoting mitochondrial homeostasis and muscle health in aging population would be to use specific activators targeting these regulators, i.e. SIRT1 agonist resveratrol and AMPK activator AICAR, which have shown beneficial effects in promoting muscle function and alleviating aging-associated muscle dysfunction in rodent models [8]. Yet, their clinical applications were hindered by uncertainties in side effects and off-targets in other organs, attesting to the urgency of uncovering novel biomarkers and potential therapeutic targets against aging-associated muscle dysfunction. Previous studies showed that miRNAs and circRNAs played important roles in the myogenic processes and muscle functions, which inspired us to investigate whether circRNAs may mediate the benefits of aerobic exercise on muscle health. Indeed, we found significant alterations in circRNAs expression pattern on the genomic scale in skeletal muscle in aging mice compared to young controls. Interestingly, GO and KEGG analysis indicated that their host genes were related to muscle contraction and development, indicating a potential cooperation between host genes and their circRNAs. Furthermore, KEGG pathway analysis of predicted target genes attributed metabolic pathway as top pathways affected by these circRNAs, unveiling a previously unappreciated role of energy metabolism in aging-associated muscle dysfunction.

Importantly, by comprehensive analysis of circRNAs profiles among young, aging and aging mice with aerobic exercise, we found that compared to young mice, circBBS9 featured decreased levels in aging mice and reversed expression after aerobic exercise and it is annotated in the circBase [22, 27]. Intriguingly, miRNA-mRNA network analysis of circBBS9 annotated top genes in metabolic pathways from KEGG analysis as its potential target genes, which supported a possible regulatory role of circBBS9 in the metabolic programs in myocytes. For example, LDHA, AMD1, DNMT3A and ACADSB were reported to regulate muscle cell proliferation and muscle atrophy [2831], while ACADS, CS, GATM, CHSY1 were related to muscle energy metabolism [3234]. Thus, aging and aerobic training, two opposite processes impact muscle functionality, may exert their impacts by modulating the circBBS9-miRNA-mRNA network. Indeed, we have validated the alternations of a few predicted circBBS9 target genes including Dnmt3a, Dad1, Gys1, Cacnale, Adcy1, Adcyap1r1 and Ctnnd1 after circBBS9 overexpression in C2C12 myotubes. Future investigations on circBBS9 in aging muscle would be informative to reveal the exact role and mechanism of its regulation, as well as the potential implication in other models of muscle dysfunction.

Kotb et al. identified and annotated many circRNAs in the skeletal muscles of Rhesus monkey via RNA-seq and verified the downregulation of mmu_circ_017332, mmu_ circ_014269, mmu_circ_015060, mmu_circ_006895, and mmu_circ_014509 during aging in monkeys [35]. In addition, previous reports have identified circFUT10, circLM07, circSVIL, circRBFOXO2, circFGFR2, circZFP609, circFGFR4 and circSNX29 were involved in myogenesis [3643]. In our study, we highlighted that circBBS9 was a potential biomarker of aging-associated muscle dysfunction and may potentially mediate the amelioration effects exerted by exercise. Of special relevance, our analysis was performed in different animal models than previous studies, and emphasized on unique changes in circRNAs in the reversal of muscle function loss with aerobic exercise, thus may result in different results.

It has to be noticed that the present study was majorly based on bioinformatic analysis of muscle samples from young and aging mice, as well as aging mice with aerobic training. Although circBBS9 was overexpressed in differentiated C2C12 myotubes and several predicted target genes, as well as genes related to mitochondrial and dystrophic functions were confirmed to be altered, suggesting the involvements of circBBS9 in muscle dysfunctions, the detailed molecular evidences of direct binding of circBBS9 with its predicted target miRNAs are warranted and in vivo application of circBBS9 in the prevention and treatment of muscle dysfunctions would be attempted in future studies. Besides, it would be interesting to examine the expression levels of circBBS9 in the circulation and its potential changes in serum/plasma under the circumstances of muscle aging and after aerobic training, which would offer circBBS9 as a novel biomarker for muscle aging.

In summary, the present study explored changes in circRNA profiles in skeletal muscle in aging and identified circBBS9 as a possible biomarker of aging-associated muscle dysfunction and potential target of aerobic exercise, which may provide new insights to the mechanism of muscle aging, and promote individualized prevention and treatment for muscle dysfunction.

Materials and Methods

Animal information and exercise training protocol

C57BL/6J mice at the age of 3-month (n=6) and 18-month (n=12, sedentary or exercise) were defined as young and aging groups, which were acquired from Shanghai Laboratory Animal Company (SLAC, Shanghai, China) and housed in pathogen-free cages at ECNU Animal Center in Shanghai. Mice were located in environment temperature at 22°C and on a 12 hours light/12 hours dark cycles with free access to water and chow diet and all of the procedures were agreed by Animal Care and Use Committee of East China Normal University (M20170316, Shanghai, China).

For exercise training, aging mice were trained with or without treadmill exercise for two months as previously reported [44] (n=6 per group). Briefly, the treadmill running assay was performed on a motorized and speed-controlled treadmill system (ZH-PT, Hangzhou, China). Mice were warmed up at a speed of 8 m/min for 5 min and then the speed was increased by steps of 2 m/min each 2 minutes until 14 m/min without inclination and sustained for 30 minutes. Mice were under aerobic training every other day for two months. The skeletal muscle tissues quadriceps femoris (Qu) were harvested from young mice, aging mice and aging mice with exercise and then frozen in liquid nitrogen and stored at -80 degrees for further studies.

Library preparation and illumina sequencing

Total RNAs were extracted from mice frozen quadriceps femoris tissues using RNAiso Plus (Takara, Japan) as products descriptions. The total RNA purity and quantity were assayed with Nano Drop 2000(USA). The RNA integrity number (RIN) ≥7 and library was constructed as manufacture’s protocols There were five microgram total RNA were used to eliminate ribosome RNA according to the description of Ribo-Zero rRNA Removal Kit (Illumina, USA). The residual RNAs were dealt with Ribonuclease R, E. coli (Epicenter, USA) to remove liner RNA and enriching circRNA, which were fragmented into oddments with RNA fragment Kit (Ambion, Austin, TX, USA). The oddments were reversely transcribed to produce cDNA. A-bases were added into ends of cDNA strand for connecting the index adaptors linked T-bases. The cDNA libraries were then stablished by PCR and pair-end sequenced (PE150) on an Illumina Hiseq 4000 machine (Genergy Bio company, Shanghai, China) according to vender’s specifications.

RNA-seq analysis and identification of CircRNAs

Cutadapt (https://cutadapt.readthedocs.io/en/stable/) was used to remove adapter sequences, primers, poly-A tails and other unwanted sequence from high-throughput sequencing reads. FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to evaluate sequence quality. circRNA analysis were exerted as described previously [45]. Bowtie2/TopHat2, TopHat-Fusion and CIRCExplorer2 softwares were used to map reads from mouse genome (Ensemble). Differential expression analysis was performed with EB-seq R package. EB-seq R package provides posterior probabilities which are adjusted for multiplicities using Benjamini-Hochberg method [46].

Gene Ontology (GO) term and KEGG pathway analysis

Gene Ontology (GO: http://www.geneontology.org/) is used to perform functional studies on gene sets [47]. Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG: https://www.kegg.jp/) database is used to understand the high-level functions and utilities of the biological system [48]. The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) is a comprehensive set of functional annotation tools for researchers to understand biological meaning behind large scale of genes [49, 50]. For analysis of cirRNA host genes and target mRNAs, GO enrichment and KEGG pathway analysis were performed using DAVID. P < 0.05 was set as the cut-off criterion.

Validation of circRNA

Total RNA was digested with RNAse R (epicenter, USA) and then reversely transcribed using random hexamers and reverse transcriptase (Takara, Japan). The quantitative PCR (qPCR) analysis was performed using divergent primers and SYBR Green (ANZY Biotech, China) on a Roche LightCycler 480 II machine. GAPDH were used to internal control for circular RNA qPCR. Besides, conventional PCR was performed using divergent and convergent primers with the templates of cDNA and gDNA respectively and PCR products were size-separated by electrophoresis with 2% of agarose gels and pictured on an UV instrument (Tiangen, China). The primers were listed in Supplementary Table 6.

Identification of sponge miRNA and miRNA-target interactions

The Miranda and RNAhybrid bioinformatics tools were utilized to predict the sponge miRNAs for circBBS9 with a maximum binding free energy of -25 kcal/mol for the sponge miRNA interaction. Next, using database Targetscan, mRNAs predicted as the targets of at least three miRNAs were overlapped and considered as mRNA targets.

Network visualization and topological analysis

We used Cytoscape software (version 3.6.1) [51] to construct and visualize the network in this study. Several topological properties such as the node degree, betweenness were analyzed using the built-in Network Analyzer tool (Department of Computational Biology and Applied Algorithmics at the Max Planck Institute for Informatics, Saarbrücken, Germany). The degree of a node is the number of edges that link to this node. Betweenness is a measure of the centrality of the node in a network, which is the number of shortest paths from each node to all others that pass through the node. The power law distribution appears as a straight line with a slope of a power exponent, which is the basis for judging the whether a random variable satisfies a power law and the power-law distribution of node degree is the main parameter used to evaluate the network topology.

Cell culture

C2C12 cells and HEK293T cells were obtained from ATCC and were cultured in DMEM containing 10% FBS, 100U/ml penicillin and 100ug/ml streptomycin. All cells were incubated at 37°C and 5% CO2 with humid atmosphere. For C2C12 cell differentiation into myotubes, culture plates were pre-incubated with 0.1% gelatin and when C2C12 cells were grown to 70% confluence, differentiated media (2% horse serum) was added and refresh every day for 5 days.

Plasmids construction and lentivirus infection

circBBS9 was inserted into a pLo-ciR vector (a generous gift from Dr. Junjie Xiao’s lab, Shanghai University). pLo-CirBBS9, psPAX2 and pMD2.G plasmids were co-transfected into HEK293T cells to generate lentiviral particles. Infectious lentivirus was harvested at 36 and 60 hours after transfection and filtered with 0.45 um filters. Differential C2C12 myotubes were infected with circBBS9 lentivirus or control lentivirus (1*109PFU) in the presence of polybrene for 48 hours. To examine circBBS9 and mRNA levels, 18S was used to internal control gene and 2-ΔΔCt methods was applied to calculated relative gene level. The primers were listed in Supplementary Table 6.

Statistical analysis

Data were analyzed and Graphics with GraphPad prism 7 version (GraphPad Software, Inc, La Jolla, CA) and R 3.3.2 version (CRNA). Data were expressed as the Mean ± SEM. Student’s t-test was used to calculate differences between two groups. p value less than 0.05 was considered as significant difference.

Author Contributions

L.X., X.G., J.L and X.M., conceived the project. M.G., H.Z, J.Y and D.W, performed biochemical and cellular experiments, F.S., Y.H. and Y.Z. established animal models and M.G., J.Q, S.W, S.X, T.H and F.S. performed RNA-seq data bioinformatic analysis. L.X. and X.M. wrote the manuscript. X.G., and J.L. edited and revised the manuscript.

Conflicts of Interest

The authors have declared no conflicts of interest.

Funding

This project is supported by funds from the National Natural Science Foundation of China (31770840 to X.M and 31800989 to L.X), the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning (TP2017042 to X.M), Shanghai Young Top Notch Talents (to X.M), Natural Science Foundation of Zhejiang Province (LY20H070003 to X.G) and ECNU public platform for Innovation (011).

References

  • 1. Christensen K, Doblhammer G, Rau R, Vaupel JW. Ageing populations: the challenges ahead. Lancet. 2009; 374:1196–208. https://doi.org/10.1016/S0140-6736(09)61460-4 [PubMed]
  • 2. Kontis V, Bennett JE, Mathers CD, Li G, Foreman K, Ezzati M. Future life expectancy in 35 industrialised countries: projections with a Bayesian model ensemble. Lancet. 2017; 389:1323–35. https://doi.org/10.1016/S0140-6736(16)32381-9 [PubMed]
  • 3. van den Beld AW, Kaufman JM, Zillikens MC, Lamberts SWJ, Egan JM, van der Lely AJ. The physiology of endocrine systems with ageing. Lancet Diabetes Endocrinol. 2018; 6:647–58. https://doi.org/10.1016/S2213-8587(18)30026-3 [PubMed]
  • 4. Boengler K, Kosiol M, Mayr M, Schulz R, Rohrbach S. Mitochondria and ageing: role in heart, skeletal muscle and adipose tissue. J Cachexia Sarcopenia Muscle. 2017; 8:349–69. https://doi.org/10.1002/jcsm.12178 [PubMed]
  • 5. Chia CW, Egan JM, Ferrucci L. Age-Related Changes in Glucose Metabolism, Hyperglycemia, and Cardiovascular Risk. Circ Res. 2018; 123:886–904. https://doi.org/10.1161/CIRCRESAHA.118.312806 [PubMed]
  • 6. Donato AJ, Machin DR, Lesniewski LA. Mechanisms of Dysfunction in the Aging Vasculature and Role in Age-Related Disease. Circ Res. 2018; 123:825–48. https://doi.org/10.1161/CIRCRESAHA.118.312563 [PubMed]
  • 7. Li L, Xiong WC, Mei L. Neuromuscular Junction Formation, Aging, and Disorders. Annu Rev Physiol. 2018; 80:159–188. https://doi.org/10.1146/annurev-physiol-022516-034255 [PubMed]
  • 8. Partridge L, Deelen J, Slagboom PE. Facing up to the global challenges of ageing. Nature. 2018; 561:45–56. https://doi.org/10.1038/s41586-018-0457-8 [PubMed]
  • 9. Fontana L, Partridge L. Promoting health and longevity through diet: from model organisms to humans. Cell. 2015; 161:106–18. https://doi.org/10.1016/j.cell.2015.02.020 [PubMed]
  • 10. Kenyon CJ. The genetics of ageing. Nature. 2010; 464:504–12. https://doi.org/10.1038/nature08980 [PubMed]
  • 11. Fontana L, Partridge L, Longo VD. Extending healthy life span—from yeast to humans. Science. 2010; 328:321–26. https://doi.org/10.1126/science.1172539 [PubMed]
  • 12. Kapahi P, Kaeberlein M, Hansen M. Dietary restriction and lifespan: lessons from invertebrate models. Ageing Res Rev. 2017; 39:3–14. https://doi.org/10.1016/j.arr.2016.12.005 [PubMed]
  • 13. Most J, Tosti V, Redman LM, Fontana L. Calorie restriction in humans: an update. Ageing Res Rev. 2017; 39:36–45. https://doi.org/10.1016/j.arr.2016.08.005 [PubMed]
  • 14. Zurlo F, Larson K, Bogardus C, Ravussin E. Skeletal muscle metabolism is a major determinant of resting energy expenditure. J Clin Invest. 1990; 86:1423–27. https://doi.org/10.1172/JCI114857 [PubMed]
  • 15. Leto D, Saltiel AR. Regulation of glucose transport by insulin: traffic control of GLUT4. Nat Rev Mol Cell Biol. 2012; 13:383–96. https://doi.org/10.1038/nrm3351 [PubMed]
  • 16. Taniguchi CM, Emanuelli B, Kahn CR. Critical nodes in signalling pathways: insights into insulin action. Nat Rev Mol Cell Biol. 2006; 7:85–96. https://doi.org/10.1038/nrm1837 [PubMed]
  • 17. Egan B, Zierath JR. Exercise metabolism and the molecular regulation of skeletal muscle adaptation. Cell Metab. 2013; 17:162–84. https://doi.org/10.1016/j.cmet.2012.12.012 [PubMed]
  • 18. Rowe GC, Safdar A, Arany Z. Running forward: new frontiers in endurance exercise biology. Circulation. 2014; 129:798–810. https://doi.org/10.1161/CIRCULATIONAHA.113.001590 [PubMed]
  • 19. Horak M, Novak J, Bienertova-Vasku J. Muscle-specific microRNAs in skeletal muscle development. Dev Biol. 2016; 410:1–13. https://doi.org/10.1016/j.ydbio.2015.12.013 [PubMed]
  • 20. Kovanda A, Režen T, Rogelj B. MicroRNA in skeletal muscle development, growth, atrophy, and disease. Wiley Interdiscip Rev RNA. 2014; 5:509–25. https://doi.org/10.1002/wrna.1227 [PubMed]
  • 21. Li J, Chan MC, Yu Y, Bei Y, Chen P, Zhou Q, Cheng L, Chen L, Ziegler O, Rowe GC, Das S, Xiao J. miR-29b contributes to multiple types of muscle atrophy. Nat Commun. 2017; 8:15201. https://doi.org/10.1038/ncomms15201 [PubMed]
  • 22. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013; 495:333–38. https://doi.org/10.1038/nature11928 [PubMed]
  • 23. Salzman J. Circular RNA Expression: Its Potential Regulation and Function. Trends Genet. 2016; 32:309–16. https://doi.org/10.1016/j.tig.2016.03.002 [PubMed]
  • 24. Han B, Chao J, Yao H. Circular RNA and its mechanisms in disease: from the bench to the clinic. Pharmacol Ther. 2018; 187:31–44. https://doi.org/10.1016/j.pharmthera.2018.01.010 [PubMed]
  • 25. Lin J, Handschin C, Spiegelman BM. Metabolic control through the PGC-1 family of transcription coactivators. Cell Metab. 2005; 1:361–70. https://doi.org/10.1016/j.cmet.2005.05.004 [PubMed]
  • 26. Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: A new star of noncoding RNAs. Cancer Lett. 2015; 365:141–48. https://doi.org/10.1016/j.canlet.2015.06.003 [PubMed]
  • 27. Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, Hanan M, Behm M, Bartok O, Ashwal-Fluss R, Herzog M, Schreyer L, Papavasileiou P, et al. Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Mol Cell. 2015; 58:870–85. https://doi.org/10.1016/j.molcel.2015.03.027 [PubMed]
  • 28. Qiu H, Xu X, Fan B, Rothschild MF, Martin Y, Liu B. Investigation of LDHA and COPB1 as candidate genes for muscle development in the MYOD1 region of pig chromosome 2. Mol Biol Rep. 2010; 37:629–36. https://doi.org/10.1007/s11033-009-9882-y [PubMed]
  • 29. Kemaladewi DU, Benjamin JS, Hyatt E, Ivakine EA, Cohn RD. Increased polyamines as protective disease modifiers in congenital muscular dystrophy. Hum Mol Genet. 2018; 27:1905–12. https://doi.org/10.1093/hmg/ddy097 [PubMed]
  • 30. Naito M, Mori M, Inagawa M, Miyata K, Hashimoto N, Tanaka S, Asahara H. Dnmt3a Regulates Proliferation of Muscle Satellite Cells via p57Kip2. PLoS Genet. 2016; 12:e1006167. https://doi.org/10.1371/journal.pgen.1006167 [PubMed]
  • 31. Chang H, Jiang S, Ma X, Peng X, Zhang J, Wang Z, Xu S, Wang H, Gao Y. Proteomic analysis reveals the distinct energy and protein metabolism characteristics involved in myofiber type conversion and resistance of atrophy in the extensor digitorum longus muscle of hibernating Daurian ground squirrels. Comp Biochem Physiol Part D Genomics Proteomics. 2018; 26:20–31. https://doi.org/10.1016/j.cbd.2018.02.002 [PubMed]
  • 32. Zoladz JA, Koziel A, Broniarek I, Woyda-Ploszczyca AM, Ogrodna K, Majerczak J, Celichowski J, Szkutnik Z, Jarmuszkiewicz W. Effect of temperature on fatty acid metabolism in skeletal muscle mitochondria of untrained and endurance-trained rats. PLoS One. 2017; 12:e0189456. https://doi.org/10.1371/journal.pone.0189456 [PubMed]
  • 33. Lee NK, MacLean HE. Polyamines, androgens, and skeletal muscle hypertrophy. J Cell Physiol. 2011; 226:1453–60. https://doi.org/10.1002/jcp.22569 [PubMed]
  • 34. Mohamed R, Dayati P, Mehr RN, Kamato D, Seif F, Babaahmadi-Rezaei H, Little PJ. Transforming growth factor-β1 mediated CHST11 and CHSY1 mRNA expression is ROS dependent in vascular smooth muscle cells. J Cell Commun Signal. 2019; 13:225–33. https://doi.org/10.1007/s12079-018-0495-x [PubMed]
  • 35. Abdelmohsen K, Panda AC, De S, Grammatikakis I, Kim J, Ding J, Noh JH, Kim KM, Mattison JA, de Cabo R, Gorospe M. Circular RNAs in monkey muscle: age-dependent changes. Aging (Albany NY). 2015; 7:903–10. https://doi.org/10.18632/aging.100834 [PubMed]
  • 36. Li H, Yang J, Wei X, Song C, Dong D, Huang Y, Lan X, Plath M, Lei C, Ma Y, Qi X, Bai Y, Chen H. CircFUT10 reduces proliferation and facilitates differentiation of myoblasts by sponging miR-133a. J Cell Physiol. 2018; 233:4643–51. https://doi.org/10.1002/jcp.26230 [PubMed]
  • 37. Wei X, Li H, Yang J, Hao D, Dong D, Huang Y, Lan X, Plath M, Lei C, Lin F, Bai Y, Chen H. Circular RNA profiling reveals an abundant circLMO7 that regulates myoblasts differentiation and survival by sponging miR-378a-3p. Cell Death Dis. 2017; 8:e3153. https://doi.org/10.1038/cddis.2017.541 [PubMed]
  • 38. Ouyang H, Chen X, Li W, Li Z, Nie Q, Zhang X. Circular RNA circSVIL Promotes Myoblast Proliferation and Differentiation by Sponging miR-203 in Chicken. Front Genet. 2018; 9:172. https://doi.org/10.3389/fgene.2018.00172 [PubMed]
  • 39. Li H, Wei X, Yang J, Dong D, Hao D, Huang Y, Lan X, Plath M, Lei C, Ma Y, Lin F, Bai Y, Chen H. circFGFR4 Promotes Differentiation of Myoblasts via Binding miR-107 to Relieve Its Inhibition of Wnt3a. Mol Ther Nucleic Acids. 2018; 11:272–83. https://doi.org/10.1016/j.omtn.2018.02.012 [PubMed]
  • 40. Peng S, Song C, Li H, Cao X, Ma Y, Wang X, Huang Y, Lan X, Lei C, Chaogetu B, Chen H. Circular RNA SNX29 Sponges miR-744 to Regulate Proliferation and Differentiation of Myoblasts by Activating the Wnt5a/Ca2+ Signaling Pathway. Mol Ther Nucleic Acids. 2019; 16:481–93. https://doi.org/10.1016/j.omtn.2019.03.009 [PubMed]
  • 41. Legnini I, Di Timoteo G, Rossi F, Morlando M, Briganti F, Sthandier O, Fatica A, Santini T, Andronache A, Wade M, Laneve P, Rajewsky N, Bozzoni I. Circ-ZNF609 Is a Circular RNA that Can Be Translated and Functions in Myogenesis. Mol Cell. 2017; 66:22–37.e9. https://doi.org/10.1016/j.molcel.2017.02.017 [PubMed]
  • 42. Wang Y, Li M, Wang Y, Liu J, Zhang M, Fang X, Chen H, Zhang C. A Zfp609 circular RNA regulates myoblast differentiation by sponging miR-194-5p. Int J Biol Macromol. 2019; 121:1308–13. https://doi.org/10.1016/j.ijbiomac.2018.09.039 [PubMed]
  • 43. Shen L, Gan M, Tang Q, Tang G, Jiang Y, Li M, Chen L, Bai L, Shuai S, Wang J, Li X, Liao K, Zhang S, Zhu L. Comprehensive Analysis of lncRNAs and circRNAs Reveals the Metabolic Specialization in Oxidative and Glycolytic Skeletal Muscles. Int J Mol Sci. 2019; 20:e2855. https://doi.org/10.3390/ijms20122855 [PubMed]
  • 44. Zhang Q, Xu L, Xia J, Wang D, Qian M, Ding S. Treatment of Diabetic Mice with a Combination of Ketogenic Diet and Aerobic Exercise via Modulations of PPARs Gene Programs. PPAR Res. 2018; 2018:4827643. https://doi.org/10.1155/2018/4827643 [PubMed]
  • 45. Lu T, Cui L, Zhou Y, Zhu C, Fan D, Gong H, Zhao Q, Zhou C, Zhao Y, Lu D, Luo J, Wang Y, Tian Q, et al. Transcriptome-wide investigation of circular RNAs in rice. RNA. 2015; 21:2076–87. https://doi.org/10.1261/rna.052282.115 [PubMed]
  • 46. Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stewart RM, Kendziorski C. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013; 29:1035–43. https://doi.org/10.1093/bioinformatics/btt087 [PubMed]
  • 47. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000; 25:25–9. https://doi.org/10.1038/75556 [PubMed]
  • 48. Altermann E, Klaenhammer TR. PathwayVoyager: pathway mapping using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. BMC Genomics. 2005; 6:60. https://doi.org/10.1186/1471-2164-6-60 [PubMed]
  • 49. Huang W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009; 37:1–13. https://doi.org/10.1093/nar/gkn923 [PubMed]
  • 50. Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009; 4:44–57. https://doi.org/10.1038/nprot.2008.211 [PubMed]
  • 51. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]