Determination of IMF content at Mth-2 and Mth-12 sheep
We measured the IMF in the longissimus dorsi of Mth-2 and Mth-12 sheep. The results show that IMF content in Mth-12 sheep (11.2 ± 0.9%) was significantly higher than in Mth-2 sheep (2.2 ± 0.006%).
Sequencing and localization of the circRNA in sheep muscle
The preliminary filtered data are listed in Additional file 1 (Table S1). The main types identified in this study were all-exon circRNA type (Mth-2: 89.36%; Mth-12: 88.56%), introns composed of lasso-type circRNA (Mth-2: 4.42%; Mth-12: 5.03 %), and circRNA in the intergenic type (Mth-2: 6.22%; Mth-12: 6.41%) (Fig. 1a and 1b). circRNA usually consists of one to four exons (Fig. 1c). It has been reported that the length of circRNA with only one exon is longer than that of circRNA with multiple exons [38]. In this study, we also found this conclusion. We list the 20 longest circRNAs in Additional file 2 (Table S2). All these are composed of one exon, most belong to the intergenic type, and are distributed on different chromosomes. A different number of circRNAs were detected on different chromosomes, with most on chromosomes 1, 2, and 3 (Fig. 1d).
Expression analysis of circRNA in muscle tissue
A total of 11,565 candidate circRNAs were identified. Mth-2 samples expressed 8,186 circRNAs, and Mth-12 samples expressed 6,683. Of these, 3,304 circRNAs were expressed in both age groups (Fig. 2a). The 30 circRNAs with the highest expression levels are listed, and the type of circRNA is marked in Table 1. Only a few belonged to ciRNA and intergenic RNAs regions. The expression of the circRNAs is shown in Fig. 2b. It can be seen from Fig. 2c that the peak gene density of circRNA expression is between 0.75 and 1.5.
Identification of differentially expressed circRNAs
Among the 104 differentially expressed circRNAs detected (Fig. 3a), 38 were upregulated, and 66 were downregulated. The overall distribution of the differentially expressed circRNAs is visually presented as a scatter plot in Fig. 3b and as a cluster heat map in Fig. 3c. The related data are listed in Additional file 3 (Table S3). Among the differentially expressed circRNAs, the source gene of circRNA3117, circRNA7648, circRNA4474, circRNA153, circRNA1817, and circRNA3096 is Titin (TTN), and the source gene of circRNA2639 is USP34. We found that the expression of circRNA829 in Mth-12 was significantly higher than in Mth-2. These circRNAs may be related to IMF, and we will study them further in the future.
Differentially expressed circRNA-hosting genes based on GO enrichment analysis
After Gene Ontology (GO) enrichment analysis, 251 groups were enriched in the biological process classification. The most significantly enriched group was regulation of transcription, DNA-templated (GO: 0006355), followed by positive regulation of GTPase activity (GO: 0043547), signal transduction (GO: 0007165), positive regulation of transcription by RNA polymerase II (GO: 0045944), and protein phosphorylation (GO: 0006468). In our research, some biological process classification GO terms related to fat metabolism were enriched, including lipid transport (GO: 0006869), negative regulation of canonical Wnt signaling pathway (GO: 0090090), lipid metabolic process (GO: 0006629), and other processes. Eighty-seven groups were significantly enriched in the cellular component classification. The most significantly enriched group was the nucleus (GO: 0005634), followed by integral component of membrane (GO: 0016021), cytosol (GO:0005829), and membrane (GO:0016020). In the molecular function classification, 114 groups were significantly enriched. The most significantly enriched group was protein binding (GO: 0005515), followed by metal ion binding (GO:0046872), ATP binding (GO:0005524), and nucleotide binding (GO:0000166). We sorted the classifications in descending order according to the number of differentially expressed genes annotated with each GO term (S gene number), and selected the top 15 in all three classification categories to draw a histogram (Fig. 4a). We also took the top 20 GO terms according to the significance of enrichment (p-value) to draw a scatter plot (Fig. 4b).
Differentially expressed circRNA-hosting genes based on the KEGG enrichment analysis
After the KEGG enrichment analysis, 95 pathways were significantly enriched. We took the top 20 KEGG terms ccording to the significance of enrichment (p-value) and drew a scatter plot (Fig. 5). More genes were enriched in Fc gamma R-mediated phagocytosis (ko04666), adherens junction (ko04520), and thyroid hormone signaling pathway (ko04919).
According to the KEGG enrichment analysis, we identified pathways with more enriched genes and pathways that might be related to IMF growth and development. These included Fc gamma R-mediated phagocytosis (ko04666), adherens junction (ko04520), thyroid hormone signaling pathway ( ko04919), fat digestion and absorption (ko04975), Notch signaling pathway (ko04330), sphingolipid metabolism (ko00600), ether lipid metabolism (ko00565), glycerolipid metabolism (ko00561), glycerophospholipid metabolism (ko00564), glutamatergic synapse (ko04724), mitogen-activated protein kinases (MAPK) signaling pathway (ko04010), and sphingolipid signaling pathway (ko04071). These 12 pathways involve ten circRNAs and their ten corresponding genes (Additional file 4: Table S4). The association between these genes and circRNAs is presented in Table 2.
Analysis of the interaction between circRNA and miRNA
CircRNA can act as a miRNA sponge, thus inhibiting miRNA from binding to target genes [39]. We selected 12 circRNAs and 18 miRNAs to predict their target relationship using TargetScan score ≥ 90 and miRanda Energy ≤ -40 (Additional file 5: Table S5), and then constructed a network diagram based on the results (Fig. 6). Among them, one circRNA could bind to multiple miRNAs, and one miRNA could bind to multiple circRNAs. We found that oar-miR-487a-5p is targeted by circRNA2748, circRNA7711, circRNA4736, and circRNA9161. We also found that circRNA7771 targets oar-miR-432_1ss23GT, oar-miR-541-5p, oar-miR-412-5p_R+2, oar-miR-1197-3p, oar-miR-323c, oar-miR-487b-5p_L-1R+1, and oar-miR-181a.
Verification of circRNA expression level
We randomly selected six differentially expressed circRNAs and detected their expression levels by quantitative real-time PCR (RT-qPCR) (Fig. 7a). The results were consistent with the trends observed in the RNA-Seq data (Fig. 7b). This shows that RNA-Seq analysis is reliable. As shown in Figure 7c, the selected circRNAs could resist digestion by RNase R, while the linear RNA (GAPDH) in the samples could not. RT-qPCR performed after RNase R digestion showed that the relative expression level of the six circRNAs did not decrease significantly. Instead, in most it has actually increased. The results show that circRNA can resist digestion by RNase R, while linear RNA cannot.