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

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

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.

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.

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.

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 5A–5C). 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).

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.

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.

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.

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*10⁹PFU) 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.