Integrated analysis of the whole-transcriptome of sheep skeletal muscle reveals the ceRNA regulation network related to muscle ber formation in sheep

Meiying Fang (  meiying@cau.edu.cn ) National Engineering Laboratory for Animal, China Agricultural University Ran Cui National Engineering Laboratory for Animal, China Agricultural University Xiaolong Kang Ningxia University Yufang Liu Hebei University of Engineering Zhen Li National Engineering Laboratory for Animal, China Agricultural University Ximing Liu National Engineering Laboratory for Animal, China Agricultural University Shuheng Chan National Engineering Laboratory for Animal, China Agricultural University Yubei Wang National Engineering Laboratory for Animal, China Agricultural University


Introduction
As a unique local sheep in Ningxia, China, Tan sheep possesses the characteristics of typical sheep in the northwestern arid and deserti ed grassland areas of China. The mutton of Tan sheep have many characteristics including tender, juicy, nutritious and unique avor. As a cross between Black headed Persian and Dorset from South Africa, Dorper sheep is famous for its rapid muscle growth for producing a desirable carcass [1]. Dutanhan sheep is a cross population of Dorper sheep, Tan sheep and small-tailed Han sheep, whose meat quality indicators are between that of Dorper sheep and Tan sheep [2].
The proportion of slow-twitch ber in sheep is positively correlated with various aspects of meat quality, such as tenderness, avor and juiciness, and as such is one of the key traits related to consumer preference [3][4][5]. Skeletal muscle is composed of two types of muscle bers that are classi ed as slow-twitch and fast-twitch ber, which show different metabolic and contractile properties. Slow-twitch bers exhibit a high oxidative capacity and a high mitochondrial content and are resistant to fatigue. In contrast, fast-twitch bers display low oxidative metabolism and low mitochondrial content. With the deepening of research, more and more regulatory factors have been identi ed to be involved in the formation of slow-twitch ber, including calcium/calmodulin-dependent protein kinase (Ca 2+ /CaMKs) signaling pathways[6-8], Ras-MAPK signaling pathway [9], Wnt signaling pathway [10], et al. However, more effective selection markers for slow-twitch ber formation participated in the sheep needs to be found. Transcriptional pro ling is a powerful approach for identifying transcripts and their expression patterns, especially more and more studies have shown that circRNAs are involved in transcriptional regulation [8] and cellular communication and signal transduction [10], and play a variety of important roles in muscle growth and development [9].
In this study, the aim of our work is to investigate the transcriptional patterns for coding genes and non-coding genes that may regulate different type of muscle ber formation in Tan sheep, Dutanhan sheep and Dorper sheep. The results would provide insight into the complex molecular mechanisms underlying transcriptional network, which contributes to increase the ratio of slow-twitch ber in actual production to improve the quality of sheep meat.

Materials And Methods
Animal Sample Collection raised in identical conditions in Ningxia, China. Sheep at the 8th month of growth were slaughtered in the same day. After slaughter, the longissimus dorsi and biceps femoris muscle samples were collected from each sheep and stored until use at −80 °C or in liquid nitrogen. Animal treatment and sample collection were approved by the Animal Welfare Committee of the State Key Laboratory of Agricultural Biotechnology, China Agricultural University, and in accordance with the Regulations on Administration of Affairs Concerning Experimental Animals, revised June 2004 (Ministry of Science and Technology).

Immuno uorescence Staining
Longissimus dorsi muscle was made into frozen sections. Tissue sections were xed for 10 min in 4% paraformaldehyde at room temperature then 0.1% TritonX-100 for 10min. After blocked for 1h using 1% BSA, tissue sections were labeled with primary antibodies (MYH7, 1/250) for overnight at 4℃. Next day, tissue sections were washed 3 times with 1xPBS and then stained with respective Alexa our 488 secondary antibody (1/250) for 1h at room temperature. Tissue sections were then washed and 10 mg/ml Hoechst was added for 15min.

RNA Extraction Library Preparation and Data Analysis
Total RNA was isolated using TRIZOL® Reagent (Invitrogen, San Diego, CA, USA) according to the manufacturer's instructions. RNA quality was assessed using 1% agarose gels. RNA purity was determined using a K5500 Spectrophotometer (Kaiao, Beijing, China). RNA integrity and concentration was assessed using the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, Foster, CA, USA).
After removed the ribosomal RNA with the Ribo-Zero rRNA kit (Epicentre, WI, USA), the RNA libraries (mRNAs, lncRNAs, and circRNAs) were generated with the mRNA-Seq Sample Preparation Kit (Illumina, CA,USA). The library quality was measured with the Agilent Bioanalyzer 2100 system (Agilent Technologies) and then sequenced using the HiSeq 4000 platform (Illumina, San Diego, CA, USA) according to the vendor's recommended protocol. Raw data (raw reads) of fastq format were obtained by removing reads containing adapter, reads on containing ploy-N and low quality reads from raw data. Reference genome (Oar_rambouillet_v1.0) and gene model annotation les were downloaded from genome website directly. Index of the reference genome was built using bowtie2 [11] and paired-end clean reads were aligned to the reference genome using HISAT2 [12]. The mapped reads of each muscle sample were assembled by StringTie [12] in a reference-based approach. After this, we evaluated the assembled transcripts using ve criteria to identify lncRNAs: (1) transcripts with exon number < 2 were removed; (2) transcripts with length ≤ 200 bp were removed; (3) known no-lncRNA annotations were removed; (4) transcripts with Fragments Per Kilobase of exon per Million fragments mapped (FPKM) < 0.5 were removed; (5) coding-non-coding-index (CNCI) [13], coding potential calculator (CPC) [14], PFAM-scan [15], phylogenetic codon substitution frequency (PhyloCSF) [16] were used to distinguish mRNAs from lncRNAs. Those transcripts that were predicted to have coding potential by all of the above tools named as candidate set of novel protein-coding transcripts, whereas those without coding potential were named as novel lncRNAs. Cuffdiff was used to calculate FPKMs of both lncRNAs and coding genes in each sample [17]. Use nd_circ[18] and CIRI2 [19] two kinds of software to lter circRNA and the results were intersected according to the position of circRNA on the chromosome. IRES nder software was used to predict whether circRNA sequences have potential IRES elements.
Refer to the standard procedure, miRNA libraries were constructed and then the quality was assessed as above cDNA libraries. The libraries were then sequenced with Illumina HiSeq™2500 system according to the vendor's recommended protocols. Clean reads were obtained after the removal of raw reads containing 5′adaptor, 3′adaptor, no insertion sequence, and poly(A) in small RNA fragments, as well as those shorter than 18nt, of known classes of RNAs (ribosomal RNAs, messenger RNAs, small nuclear RNAs, transfer RNAs, small nucleolar RNA, small cytoplasmic RNAs, and repeats), and of low quality. Bedtools (https://bedtools.readthedocs.io/) was used to search for known miRNAs by matching them to entries in miRBase20.0 (http://www.mirbase.org/). After excluding reads that were mapped to known miRNAs, miRDeep2 [20] was used to analyze the remaining reads to predict novel miRNAs. The expression amount of known and new miRNA and circRNA in each sample was counted, and the expression quantity was normalized by TPM [21].

GO and KEGG enrichment analysis
DESeq [22]was used for differential gene expression analysis of transcripts among three groups. only the transcripts with signi cant differences (P-values < 0.05) and | log 2 FC| > 1 were classi ed as DE transcripts. Gene ontology (GO) enrichment analysis of DE mRNAs, DE circRNAs parental genes , DE lncRNAs and DE miRNAs target genes were conducted using the GO-seq R package [23], correcting for gene length bias. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system [24], such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. KOBAS v3.0 [25] was used to test the statistical signi cance for the enrichment of DEGs or DE circRNAs parental genes or targets of DE lncRNA and DE miRNAs in KEGG pathways. Only the p< 0.05 was considered statistically signi cant.
LncRNA/circRNA-miRNA-mRNA Network Analysis For the cis-acting prediction, coding genes 10kb upstream and downstream of lncRNA were searched as cis-acting regulatory targets. For the trans-acting prediction, the expression of the lncRNA was determined to be not related with the location of the mRNA but co-expressed with it.
The expressed correlation was calculated between lncRNAs and coding genes using Pearson's correlation coe cients (r > 0.90 or r < −0.90) as a classi er. The prediction of miRNA target genes was performed using miRanda [26], PITA [27] and RNAhybrid [28]. CeRNA hypothesis RNA transcripts can crosstalk by competing for common miRNAs, with miRNA response elements as the foundation of this interaction [29].
lncRNA/circRNA-miRNA-mRNA pairs were constructed with lncRNA/circRNA as decoy, miRNA as the core, and mRNA as the target. The lncRNA, circRNA, miRNA and mRNA interactions were constructed and visualized using Cytoscape [30].

Quantitative Polymerase Chain Reaction
To evaluate the reliability of the transcript expression data obtained by RNA-Seq. Quantitative Polymerase Chain Reaction (qPCR) were carried out with the longissimus dorsi tissue. Using FastKing gDNA Dispelling RT SuperMix and SuperReal PreMix Plus (Tiangen, China), experiments were conducted for mRNA detection. Using miRcute Plus miRNA rst-strand cDNA Kit and miRcute Plus miRNA qPCR Kit (Tiangen, China), experiments were conducted for miRNA detection. Using RT primers and Custom gene qRT-PCR Quantitation Kit (GenePharma, China), experiments were conducted for circRNA detection. To verify the back-splicing junction of circRNA, the RNase R-and control reaction (without RNase R) were prepared and reverse transcribed into cDNA to amplify each circRNA with primers following previously method [31,32], then the products were sequenced to nd the back-splicing sites. The GAPDH and U6 small nuclear RNA genes were selected as the endogenous control genes (all primers are shown in Table S12). All qPCR validations were carried out with three biological replicates and triplicate reactions for each sample. After ampli cation, the products were con rmed by agarose gel electrophoresis and Sanger sequencing, the relative transcript abundance was calculated using 2 -ΔΔCt method.

Statistical Analysis
Data were expressed as means ± standard deviation (SD). Signi cance was analyzed using one-way analysis of variance (ANOVA) to test homogeneity of variances via Levene's test, followed by Student's t-test. Calculations were conducted using SAS version 9.0 (SAS, Cary, NC, USA). Differences were considered to be statistically signi cant for P-values < 0.05.

Results
The ratio of fast and slow muscles ber in longissimus dorsi tissue in different populations The results of immuno uorescence staining showed that the proportion of slow-twitch bers in longissimus dorsi tissue in Tan sheep was about 11.247%, and that in Dorper sheep was about 5.408%. Compared with Dorper sheep, the proportion of slow-twitch bers was signi cantly higher in Tan sheep(P<0.05), and the increase ratio accounted for about 51.9% of Tan sheep (Fig.1).

Overview of RNA sequencing
To assess the genes involved in muscle ber formation, longissimus dorsi (LD) and biceps femoris (BF) tissues were collected from the 8month-old Dorper sheep, Tan sheep and Dutanhan sheep for the whole-transcriptome pro ling of all mRNAs and noncoding RNAs. For the RNA sequencing libraries, an average of 95.74 million clean reads were obtained from the 18 samples tested, and 83.55-90.07% of these reads were uniquely aligned to the reference genome Oar_rambouillet_v1.0. All 18 samples had at least 85% reads equal to or exceeding Q30 (Table   S1-1). In addition, for the small RNA-Seq libraries, an average of 16.02 million clean reads were obtained (Table S1-2). Moreover, 107 mRNAs were differently expressed, which included 47 BF tissue-speci c and 37 LD-speci c mRNAs, respectively. Apart from mRNAs, 224 known lncRNAs and 22373 novel lncRNAs were identi ed from the RNA-seq data by blasting to lncRNAs in the NCBI database and performing CNCI, CPC2 and PfamScan analysis (Fig.S1). Additionally, 737 lncRNAs were differently expressed, of which 216 BF tissue-speci c and 242 LDspeci c. Of these DE lncRNAs, 37.9% were expressed in both types of tissues. Furthermore, 238 miRNAs were identi ed and 32 miRNAs were differently expressed, of which 22 were BF tissue-speci c and 4 were LD-speci c. Finally, 3648 circRNAs were identi ed and 348 circRNAs were differently expressed, of which 188 and 135 were found to be BF tissue-speci c and LD-speci c, respectively (

DE mRNAs related to muscle bers formation
First, DEGs between the LD and BF tissues from 8-month-old Dorper sheep, Tan sheep and Dutanhan sheep were screened, where upon 149 DEGs were found, and their log 2 FC values were presented as Volcano Plot pictures, which ranged from −20 to 20. Among all of these mRNAs, 25 (13 up-regulated, 12 down-regulated) and 35 (24 up-regulated, 11 down-regulated) DE mRNAs were identi ed in the D_L/T_L and D_B/T_B, 32 and 42 were identi ed in DTH_L/T_L and DTH_L/D_L, and 21 and 41 were identi ed in DTH_B/T_B and DTH_B/D_B, respectively (Fig.2a). The amount of DE mRNAs in BF and LD tissues was basically the same, implying that they might participate in the equally important response process under meat quality.
To explore the functions of these DE mRNAs, GO annotation and KEGG enrichment analysis were performed (Fig.2b and c). GO annotation analysis showed that these genes were involved in calcium ion binding (GO:0005509), cellular response to calcium ion (GO:0071277), ATPase inhibitor activity (GO:0042030), oxidation-reduction process (GO:0055114). KEGG enrichment analysis showed that muscle ber formation, metabolic pathways, cAMP signaling pathway were signi cantly enriched, and these terms of DE mRNAs possible have primarily participated in meat quality. Among them, DEGs included protein kinase AMP-activated catalytic subunit alpha 2 (PRKAA2, log 2 Fold change (log 2 FC=8.11), EF-hand domain family member B (EFHB, log 2 FC=2.43), L-lactate dehydrogenase A chain (LDHA, log 2 FC=3.35), folliculin interacting protein 2 (FNIP2, log 2 FC=3.05), which are related to pathways associated with muscle ber formation: the cellular response to calcium ion, calcium ion binding, oxidation-reduction process and ATPase inhibitor activity, respectively (Table S4).
In addition, 3 DEGs(ZIC1 MYL3 ASAH1) were co-differentially expressed in both LD and BF tissues (Fig.S3a), all of which had known functions in muscle bers growth and development. These differentially expressed transcripts may be related to the different function during the development of these two tissues, which are worthy of future analysis.

DE lncRNAs related to muscle bers formation
Comparison of the genomic characterizations of the lncRNAs with mRNAs showed that their transcripts were similar in length distribution; for exon number, a higher number of lncRNAs had 2 to 4 exons; in addition, lncRNAs had a shorter ORF length and lower FPKM value ( Fig.3a  To reveal the potential functions of the 838 identi ed lncRNAs in muscle ber formation, two independent algorithms-cis (genomic location) and trans (expression correlation) -were performed to predict the target genes of the lncRNAs. In total, 720 were targets of 2154 cis-  .3 d and e). As a preliminary exploration of the functional implications of the DELs across genomes, we investigated whether lncRNAs were co-regulated with the DEGs during muscle ber formation. Interestingly, in both LD and BF tissues, we observed that the cis-lncRNA LNC_002028 targeted EFHB as its DE target gene, whereas the three trans-lncRNAs LNC_003248, LNC_003249, LNC_013114 targeted PRKAA2 as their DE cotarget genes. In addition, 4 DELs(LNC_003158 LNC_015291 LNC_015292 LNC_015293) were codifferentially expressed in both LD and BF tissues (Fig.S3b).

DE circRNAs related to muscle bers formation
A total of 3651 circRNAs were identi ed in this study, and circRNAs were de ned with at least two independent reads crossing the backsplice junction.Among the 3651 circRNAs we studied, circRNAs were more likely to be produced from exons, 3141 circRNAs were generated from exons of 1810 genes, 197 circRNAs were generated from intronic regions of 175 genes, and the rest were intergenic circRNA according to our statistics (Fig. 5a).CircRNAs were mostly located on chromosome 1 and ranged in length from 192nt to 93278nt (Fig. 5b). Read counts across backsplice junction sites were normalized to transcripts per million (TPM) (Fig. 5c). The source sites of circRNAs (including exons, introns, and intergenic regions) were different, but no signi cant changes were found in the lengths of circRNAs from different sources (Fig. 5d).In addition, the expression correlation between circRNAs and their host genes was also detected based on the linear splicing sites of the parental gene and the reverse splicing sites of the circRNAs. The Pearson correlation coe cient between circRNAs and their host genes was 0.046, indicating a very low correlation (Fig. 5e).Then, the coe cient of variance of each circRNA and its parental gene was calculated as a measure of the variance of the expression level, and it was found that circRNAs may be easily regulated by multiple factors and therefore have higher expression variability (Fig. 5f).Interestingly, it was found that circRNAs had extremely high potential to encode proteins, and 4611 internal ribosome entry sites (IRES) were predicted, and the highest score could be as high as 0.997 (Table S10).  5 g and i, TableS11). In addition, novel_circ_0009997 were codifferentially expressed in both LD and BF tissues (Fig.S3d).

Construction of the ceRNA regulatory network
It has been shown that mRNAs, lncRNAs, and circRNAs may act as ceRNAs, which regulate gene function via miRNA in various processes, suggesting that ceRNAs and their miRNAs may be co-regulated in muscle ber formation. On the basis of the data of the DE mRNA, lncRNA, circRNA, and miRNA transcripts, a ceRNA network was constructed consisting of 4 DEGs, 21 DEMs, 115 DECS, and 26 DELs, and obtained 70 mRNA-miRNA, 31 lncRNA-miRNA, and 158 circRNA-miRNA interaction relationships (Fig.6, Fig.S2). Within the network, it was found that PRKAA2, EFHB, FNIP2 and LDHA may be the crucial genes mediated by noncoding RNAs for regulating muscle ber formation. The functional prediction results showed that many of these DE circRNAs, such as novel_circ_0017336, novel_circ_0007039, novel_circ_0017132, competing to bind with miR-23a, might act as ceRNAs to control PRKAA2, EFHB and FNIP2 level, respectively.

Validation of gene expression by RT-qPCR
Validation of the RNA-seq results was carried out using the qPCR for 4 DEGS (PRKAA2, EGR1, MYL3 and FBXL5), 5 DEMs (miR-26a, miR-308-3p, miR-409-3p, miR-495-3p and miR-23a), and 4 DECs (circ_0017336, circ_0017132, circ_0017430, and circ_0007039) (Fig.7 a, b and d). The expression of these selected transcripts was signi cantly different in the LD tissue in Tan sheep and Dorper sheep, with the expression patterns being highly consistent with those obtained by RNA-Seq. It's worth mentioning that the backsplice junctions of circRNAs were con rmed before the validation, due to the circular structure, the circRNAs were more resistant to digestion by RNase R treatment, and the back-splicing sites were veri ed by Sanger sequencing (Fig.7c). The results indicated the high reproducibility and reliability of the gene expression pro les obtained in our study.

Discussion
As is widely known, the muscle ber formation is one of the polygenic traits in sheep that is an important determinant of meat quality characteristics. In general, the diameter of slow-twitch ber is thinner than that of fast-twitch ber, while the density of slow-twitch ber is larger than that of fast-twitch ber, and its proportion is positively correlated with meat quality. There are some differences in histological characteristics of muscle bers among different populations of sheep. In previous studies, compared LD tissues with Dutanhan sheep and Dorper sheep, Tan sheep had the thinnest muscle ber diameter, while the muscle ber density of Tan sheep was signi cantly greater than that of Dorper sheep and Dutanhan sheep [2]. In this study, we calculated the ratio of fast and slow muscle bers in LD tissue by immuno uorescence staining, whereupon it was found that Tan sheep has a higher proportion of slow muscle bers than Dorper sheep. This prompted us to choose LD and BF two muscle tissues in different groups for the whole transcriptome pro le analysis to identi ed DEGs and DE non-coding RNAs. These changes may be closely related to the formation of slow muscle bers in sheep.
The sheep used in this study are 8-month-old population that have attained consistency in appearance, reproductive and production performances. 139 DEGs were identi ed in the muscle tissues of three groups, many of which have known functions in muscle ber formation. For example, PRKAA2 encodes the AMPKa2 subunit, and the AMPKa2 subunit is a catalytic subunit of AMPK, which participates in the activation of AMPK protein. AMPK family gene polymorphism plays a role in skeletal muscle gene expression differences and effects on meat quality [7,[33][34][35]. The expression level of PRKAA2 in Tan sheep is higher than that in Dorper sheep, which may also promote the formation of slow muscle ber by increasing the level of energy metabolism. Moreover, the expression of FNIP2 was in consistent with the increased trend of slow muscle ber, whereas EFHB, LDHA showed the opposite trend, indicating that the former gene acts as positive regulators in slow muscle ber formation whereas the latter genes act as negative regulator in this process. FNIP2 is involved in the regulation of ATPase inhibitor activity. It has been reported that FNIP may regulate AMPK pathway via direct or indirect pathways, thereby regulating mitochondrial biosynthesis and muscle ber types. EFHB participates in the regulation of calcium ion binding, the transcription factors of mitochondrial nuclear-encoded genes are activated in skeletal muscle by Ca 2+ -dependent transport pathways, and then cause the formation of muscle bers. LDHA is involved in the regulation of oxidation-reduction process and lactate dehydrogenase activity. Besides, 3 DEGs (ZIC1, MYL3, ASAH1), 4 DELs (LNC_003158, LNC_015291, LNC_015292, LNC_015293), 1 DEMs(oar-miR-10b) and 1 DECs(novel_circ_0009997) were co-differentially expressed in both LD and BF tissues, many of which had known functions in muscle bers growth and development. For instance, studies have shown that ZIC1 is differently expressed in longissimus dorsi and semitendinosus of the Wuzhishan pig, and it can play a role in Myf5 regulation and somite myogenesis[36]. MYL3, myosin light chain 3, has been identi ed as candidate function gene in muscle development in Dorper and small-tailed Han sheep [37]. While mutations in ASAH1 will cause the occurrence of muscular atrophy in rare diseases [38]. These differential expression transcripts may be related to the different function of these genes during the development of these two tissues, which are also worthy of future analysis.
Although biological processes, cellular composition, and molecular functions nally depend on the expression of protein-coding genes, over the last twenty years, the concept has been transformed that the expression of protein-coding genes in turn depend on an intricate transcriptional network, including a variety of non-coding RNAs [39][40][41]. A very abundant classe of ncRNAs (~270,000 transcripts) are lncRNAs, which are characterized by transcripts of 200nt or longer. They share many features with mRNA, including polyadenylation, 5 capping and splicing, but most are not translated into protein. lncRNAs could play a major role in many physiological and pathological conditions through trans and cis activities, and prediction of lncRNA target genes by these two independent algorithms as well as bioinformatics analysis will help to identify which term lncRNAs are participated in and further showing potential function. Our results showed that lncRNAs obtained in this study may regulate pathways related to muscle ber formation via trans and cis activities (such as skeletal muscle ber development, Negative regulation of glycolysis, ATPase activity). In addition, our results also indicated that the cis-lncRNA, LNC_002028, and three trans-lncRNAs, LNC_003248, LNC_003249, and LNC_013114, also regulated two important candidate genes, EFHB and PRKAA2, respectively. On the basis of these results, we suspected that one of the main roles of these 4 lncRNAs is to regulate the muscle ber formation in sheep, and further highlighting of their detailed mechanisms in this process would be fertile ground for future investigation.
Besides, miRNAs are small (~22nt) trans-acting functional RNAs involved in silencing and post-transcriptional regulation of gene expression [39]. The results on the DEMs and DECs in muscle tissue indicated that they were signi cantly enriched in muscle ber formation pathways, for instance, FOXO signaling pathway, Ca 2+ signaling pathway, AMPK signaling pathway, contraction ber, glycogen metabolic process. From the constructed ceRNA network, we found several miRNAs, including substantially down-regulated miR-409-3p and up-regulated miR-23a, miR-26a, miR-27a, miR-382-3p, novel-67 and miR-495-3p around PRKAA2. Studies have shown that miR-23a down-regulates the translation of mabbx/atrogin-1 and MuRF1 (mediated atrophic protein degradation) to promote muscle atrophy [40]; meanwhile, miR-23a inhibits myocyte differentiation by down-regulating MYHC gene expression[41]; A similar research of miR-23a was carried in Angus cattle at the 80 fetus (90 days old) and adult (24 months old) stages, circMYBPC1 was identi ed differently expressed and interacted with muscle genes MYHC through directly binding miR-23a[42]; miR-26a promotes myoblast differentiation by down-regulating the expression of its target gene Ezh2[43]. Our miRNA-mRNA interaction analyses identi ed a much large number of predicted miRNA-mRNA pairs in the comparison of Dorper sheep, Dutanhan sheep and Tan sheep, suggesting that Tan sheep have a greater ability to regulate the slow muscle ber forming environment. CircRNA is a novel type of non-coding RNA, which exhibits high stability, abundance, tissue/stage speci city, and evolutionary conservation [9]. All of these 3651 DECs are new circRNAs that have not been discovered before in sheep, and there is little information or data obtained from published databases. In our study, the functional prediction results also showed 11 DE circRNAs, including novel_circ_0017336, novel_circ_0007039, novel_circ_0017132, competing to bind with miR-23a, might act as ceRNAs to control PRKAA2, EFHB and FNIP2 level, respectively, which side veri es the importance of miR-23a not only in its own function, but also acting as a bridge in the ceRNA mechanism of muscle ber formation. GO and KEGG pathway analysis were constructed to the target genes involved in circRNA-or lncRNA-associated ceRNA. Interestingly, GO annotation showed that associated genes are all enriched in the same term: regulation of skeletal muscle tissue growth. Further studies on ANKRD, UCP1, ADRB3 and other genes in this pathway are also necessary.
In this study, strict constraints were applied to select the most possible circRNA-or lncRNA-associated ceRNA networks that are involved in muscle bers formation. For example, we selected ceRNA networks that participated in the control of crucial muscle gene related to signaling pathway. 4 DEGs, 115 DECs, and 26 DELs were connected with one another through the 21 DEMs, a total of 70 mRNA-miRNA, 31 lncRNA-miRNA, and 158 circRNA-miRNA interaction relationships were obtained. This indicated that muscle ber formation in sheep is jointly regulated by these factors and results from a balance level of transcripts expression. We speculate that most of these factors may play a key role in the regulation of muscle ber formation, these possibilities are worthy of further research efforts.

Conclusion
Our review gives a complete scene of the distinctions in the entire transcriptome pro les of LD and BF tissues among three gatherings in sheep. 4 DEGs were recognized connected with muscle ber arrangement among various gatherings, including PRKAA2, EFHB, FNIP2, LDHA.
Moreover, a few DE lncRNAs, miRNAs, and circRNAs were found in muscle tissue rmly connected with a few muscle ber arrangement pathways. Based on the differently expressed records, a ceRNA regulatory network was developed which contained 4 DEGS, 115 DECS, 26 DELS, 21 DEMS and 259 connections. Our discoveries have distinguished potential regulators and molecular regulatory networks that might be associated with muscle ber arrangement in sheep, and gave an establishment to future examinations on the molecular mechanisms underlying muscle ber development.   Figure 1 The difference between the slow-twitch ber ratio of the longissimus dorsi muscle of Tan sheep and Dorper sheep was observed and counted by frozen section MyH7 immuno uorescence staining.     LncRNA/circRNA-miRNA-mRNA interaction network. ovals, triangle, parallelogram and V type represent mRNAs, miRNAs, lncRNAs and circRNAs, respectively.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.