Overview of circRNA profiles in LD and Sol skeletal muscle
To determine the identity, abundance and differential expression of circRNAs in different types of skeletal fibers, we profiled circRNAs in slow-type-enriched Sol and fast-type-enriched LD of large white pigs. Ribo-depleted total RNA-Seq libraries were prepared and sequenced according to the flow chart (Fig. 1A). We identified a considerable number of RNAs in LD and Sol muscles when considering total RNA libraries (Table 1). A total of 40757 candidate circRNAs were identified by at least one unique back-spliced read. We found that 15627 and 14742 circRNAs were specifically expressed in LD and Sol muscles respectively, and 10388 were expressed in both muscles (Fig. 1B). For both the LD and Sol samples, the circRNAs are mostly comprised of exonic, intronic and intergenic sequences, as determined based on their mapping to the genome, whereas a smaller fraction of circRNAs also contains upstream, downstream and UTRs (Fig. 1C).
Identification and characteristics of circRNAs expressed in LD and Sol
To further estimate the characteristics of circRNAs, we analyzed the circRNA sequences. Standard metrics to characterize the length of circRNAs detected in this study (minimum, maximum, mean, median, and total length) are provided in Table 2. Most circRNAs were no more than 1000 nucleotides (nt), and the median length was 439 nt (Fig. 2A). According to their host gene location, the 40757 circRNAs were widely distributed across chromosomes except for the Y chromosome, and the number of circRNAs increased with chromosome length (Fig. 2B). According to the histogram depicting flanking intron lengths, the length of flanking intron regions of most circRNAs was approximately 103-104 nt, indicating that long flanking introns are necessary for circRNA formation (Fig. 2C-D). A total of 21193 (52.0%) circRNAs were composed of one to eight exons, among which 4606 (11.3%) contained one exon. Furthermore, 12204 (29.9%) circRNAs contained more than fifteen exons (Fig. 2E). We identified 40757 circRNAs generated from 8009 host genes, indicating that most genes can produce multiple circRNAs. There were 1007 (12.6%) circRNA-producing genes that generated a single circRNA, whereas 26.6% produced more than fifteen (Fig. 2F).
Analysis of differentially expressed circRNA between LD and Sol
Based on the analysis of circRNA expression, we found 181 significantly differentially expressed circRNAs when comparing the libraries derived from the LD and Sol (Additional File). To further explore the potential functions of circRNA, we constructed a clustered heatmap (Fig. 3A). Although several circRNAs showed the same expression pattern, most of them were differentially expressed between the LD and Sol. Based on the expression levels of circRNAs in paired samples (LD to Sol ratio), 90 circRNAs were downregulated, and 91 circRNAs were upregulated in the LD (Fig. 3B).
Construction of circRNA-miRNA interaction network
A previous study has shown that circRNAs composed of only exons can play a role in the cytoplasm by sponging microRNA (miRNA). We randomly selected 11 differentially expressed exonic circRNAs and predicted the miRNAs they might sponge by miRanda, RNAhybrid and TargetScan. Then, we constructed an interactive network between these 11 circRNAs and the sponged miRNAs (Fig. 4). The whole interaction network included 11 circRNAs, 111 miRNAs and 146 relationship lines.
KEGG pathway enrichment analysis
To further analyze the regulatory mechanisms of differentially expressed circRNAs in different types of skeletal muscle fibers, we predicted the sponge miRNAs of the differentially expressed circRNAs and analyzed the target genes of miRNAs. Then, we performed pathway enrichment analysis of the predicted target genes (Fig. 5). A total of 30 signaling pathways were enriched for these genes, among which metabolic pathways were the most abundant pathway. Markedly, some pathways were closely related to muscle fiber development, resulting in the transformation of muscle fiber types, such as the AMPK signaling pathway, FoxO signaling pathway, and PI3K-Akt signaling pathway.
Validation of differentially expressed circRNAs by RT-qPCR
To validate the data reliability of differentially expressed circRNAs detected by sequencing, we randomly selected 9 differentially expressed circRNAs and amplified their junction regions using specific RT-qPCR primers. RT-qPCR revealed two upregulated and seven downregulated circRNAs, as shown in Fig. 6, which clearly showed that expression patterns of the nine selected circRNAs were consistent with the results of RNA-Seq. These results demonstrated the accuracy and reliability of high-throughput sequencing.
Overexpression of circMYLK4 promotes the development of slow muscle fibers
We analyzed novel circRNAs in many aspects and selected circMYLK4, which was named after its host gene MYLK4 located on chromosome 7, for its high expression in Sol for further analysis. To investigate whether this circRNA could be a regulatory factor during muscle fiber development, we injected piglets with circMYLK4-AAV. Then, we collected gastrocnemius (GAS), semimembranosus (SB), semitendinosus (SM) and soleus (Sol) samples for analysis 4 weeks after injection. The expression level of circMYLK4 in AAV group was significantly higher than that in control group (Fig. 7A). The results of the meat color assessment showed that circMYLK4 reduced the brightness of meat and improved the redness and yellowness of meat (Fig. 7B). In support of this phenomenon, circMYLK4 overexpression significantly increased expression levels of the slow muscle marker gene MyHC I, which was detectible at both the mRNA and protein levels (Fig. 7C-D). In addition, we also measured the mRNA of fiber type-related genes (Tnnt1, Tnni1, and Tnnc1), mitochondrial biogenesis genes (PGC-1α and Tfam), and mitochondrial respiratory chain genes (Atp5o, Cycs, and Ndufa5) by RT-qPCR. As expected, circMYLK4 overexpression increased the expression of the above genes (Fig. 7E-H). Collectively, these results established that circMYLK4 promoted the development of slow muscle fibers by regulating the expression of mitochondrial-related factors.