Integrated analysis of the whole transcriptome of skeletal muscle reveals the ceRNA regulatory network related to the formation of muscle fibers in Tan sheep

Meat quality is highly influenced by the kind of muscle fiber, and it can be significantly improved by increasing the percentage of slow-twitch fibers. It is still not known which genes control the formation of muscle fibers or how those genes control the process of forming in sheep until now. In this study, we used high-throughput RNA sequencing to assess the expression profiles of coding and noncoding RNAs in muscle tissue of Tan sheep and Dorper sheep. To investigate the molecular processes involved in the formation of muscle fibers, we collected two different muscle tissues, longissimus dorsi and biceps femoris, from Tan sheep and Dorper sheep. The longissimus dorsi of Tan sheep and Dorper sheep displayed significantly differential expression levels for 214 lncRNAs, 25 mRNAs, 4 miRNAs, and 91 circRNAs. Similarly, 172 lncRNAs, 35 mRNAs, 12 miRNAs, and 95 circRNAs were differentially expressed in the biceps femoris of Tan sheep and Dorper sheep according to the expression profiling. GO and KEGG annotation revealed that these differentially expressed genes and noncoding RNAs were related to pathways of the formation of muscle fiber, such as the Ca²⁺, FoxO, and AMPK signaling pathways. Several key genes are involved in the formation of muscle fibers, including ACACB, ATP6V0A1, ASAH1, EFHB, MYL3, C1QTNF7, SFSWAP, and FBXL5. RT-qPCR verified that the expression patterns of randomly selected differentially expressed transcripts were highly consistent with those obtained by RNA sequencing. A total of 10 lncRNAs, 12 miRNAs, 20 circRNAs, and 19 genes formed lncRNA/circRNA-miRNA-gene networks, indicating that the formation of muscle fiber in Tan sheep is controlled by intricate regulatory networks of coding and noncoding genes. Our findings suggested that specific ceRNA subnetworks, such as circ_0017336-miR-23a-FBXL5, may be critical in the regulation of the development of muscle fibers, offering a valuable resource for future study of the development of muscle fibers in this animal species. The findings increase our understanding of the variety in how muscle fibers originate in various domestic animals and lay the groundwork for future research into new systems that regulate the development of muscle.

Introduction

Meat quality is an important economic trait, and the ratio of muscle fibers is one of the most important factors affecting meat quality traits, which can be significantly improved by increasing the proportion of slow-twitch fibers. Skeletal muscle is composed of two types of muscle fibers that are classified as slow-twitch and fast-twitch fibers, which have different metabolic and contractile properties. Compared with fast-twitch fibers, slow-twitch fibers have higher oxidative metabolism capacity and higher mitochondrial content, so their proportions are positively correlated with various aspects of meat quality, such as tenderness, flavor, and juiciness (Doerr, 1971; RyuandB and Kim, 2005; Su et al., 2013). Ningxia Tan sheep are an important breed for fur in China (Kang et al., 2013; Liu et al., 2018), but with the decline of the fur market, they have gradually been bred for meat. Previous studies have shown that some meat quality indicators of Tan sheep were measured, such as pH, hydraulic power, meat color, shear force, cooked meat percentage, muscle fiber density, and diameter, while the relationship between the proportion of muscle fibers and meat quality traits in Tan sheep requires a more precise interpretation; in this study, we compared the proportions of slow-twitch fibers in Tan sheep and Dorper sheep. Studies have also found the heritability of the muscle fiber content ranged from 0.20 to 0.68 (Larzul et al., 1999), and molecular breeding is a key method used for the improvement of slow-twitch fiber content. However, no effective molecular markers for slow-twitch fiber content selection practices in sheep have yet been found.

A growing number of studies have revealed that noncoding RNAs are involved in transcriptional regulation (Kotb et al., 2015), cellular communication, and signal transduction (Kristensen et al., 2019) and play a variety of significant roles in muscle growth and development (Hong et al., 2019). Competing endogenous RNA (ceRNA) is a common physiological mechanism by which lncRNAs and circRNAs compete with miRNAs for binding (Salmena et al., 2011). Previous studies have demonstrated noncoding RNA and ceRNA regulation of myoblast differentiation and the formation of muscle fibers in mouse, chicken, and bovine species. For example, the intronic lncRNA SYISL recruits polycomb repressive complex 2 in C2C12 cells to suppress the expression of the fast-twitch fiber marker gene and myoblast differentiation (Jin et al., 2018). CircPTPN4 can function as a ceRNA to regulate nicotinamide phosphoribosyltransferase (NAMPT) expression by sponging miR-499-3p, thus participating in muscle fiber type switching in chicken primary myoblasts (Cai et al., 2022). The longissimus dorsi in cattle was shown to express CircMYBPC1 differently during adult and embryonic phases of development, and CircMYBPC1 interacted with the muscle gene MYHC by directly binding to miR-23a (Chen et al., 2021). Although many studies related to the formation of muscle fibers have been published, the detailed molecular mechanism remains unclear, especially in noncoding RNA and ceRNA in sheep.

The molecular process of the formation of muscle fibers in Tan sheep has not yet undergone a comprehensive examination. Therefore, clarifying the roles of the DEGs and noncoding RNAs can contribute to understanding the changes in the whole transcriptomes associated with the formation of muscle fibers in different sheep groups. In this study, the muscle fiber contents of longissimus dorsi and biceps femoris in Tan sheep and Dorper sheep were determined, and whole-transcriptome profiles in different groups (Tan sheep and Dorper sheep) and different tissues (longissimus dorsi and biceps femoris) of the same groups were obtained to compare the DEGs and noncoding RNAs. Finally, we built a thorough competing endogenous RNA (ceRNA) network to identify the genes that are most likely to be involved in the formation of muscle fibers. In-depth knowledge of the molecular processes underlying the development of muscle fibers in Chinese Tan sheep is provided in this report.

Materials and methods

Animal sample collection

Animal care and sample collection were authorized by the State Key Laboratory of Agricultural Biotechnology of the Agricultural University of China’s Animal Welfare Committee and carried out by the June 2004 amendments to the Regulations on the Control of Experimental Animals (Ministry of Science and Technology). Three unrelated male Tan sheep were used in the experimental group, while three unrelated Dorper sheep were used in the control group. Both groups of sheep were bred in Ningxia, China, under identical conditions until they were 8 months old, at which point they were killed on the same day. From these six sheep, samples of the longissimus dorsi and biceps femoris were taken and kept in liquid nitrogen or at −80°C until needed.

Immunofluorescence staining experiment

To examine whether the proportion of slow-twitch fibers differed among different sheep groups, we performed immunofluorescence using a slow-twitch marker protein (MYH7 antibody). The fluorescein-labeled MYH7 antibody was used as a probe to detect the target antigen in the muscle tissue and to identify the location of the antigen. First, the frozen tissue of the longissimus dorsi and biceps femoris was sectioned. Slices of tissue were first fixed for 20 min at room temperature in 4% paraformaldehyde, followed by 20 min of penetration in 0.1% Triton X-100. Tissue sections were tagged overnight with an antibody (MYH7, 1/500) at 4°C following an hour of closure with 1% BSA. Tissue slices were washed three times with ×1 PBS the next day and then stained with Alexa Fluor 488-conjugated antibody (1/1,000) at room temperature for 1 h. Then, the tissue slices were washed with ×1 PBS, and 10 mg/ml Hoechst was added for 10 min.

RNA extraction and sequencing

Following the manufacturer’s instructions, total RNA was extracted using TRIzol Reagent (Invitrogen, United States). RNA quality was evaluated using agarose gels at a 1% concentration. A K5500 spectrophotometer was used to measure the purity of the RNA (Kaiao, China). The Bioanalyzer 2100 system and the RNA Nano 6000 Assay Kit were used to measure RNA integrity and concentration (Agilent Technologies, United States). The mRNA/lncRNA/circRNA library was built using 5 µg of RNA from each sample. The HiSeq 4000 sequencing platform (Novogene, China) was used to build mRNA/lncRNA/circRNA sequencing libraries, and 150 bp paired-end reads were produced by the manufacturer’s instructions.

To create the small RNA library, 3 μg of the sample’s total RNA were employed. Using the NEBNext^® Multiplex Small RNA Library Prep Set for Illumina^® (NEB, United States) by the manufacturer’s instructions, sequencing libraries were built. Using DNA High Sensitivity Chips, the Agilent Bioanalyzer 2100 system assessed the sequencing library’s quality. The HiSeq 2500 sequencing platform (Novogene, China) was used to build small RNA sequencing libraries, and 50 bp single-end reads were produced. Fastq formatted raw data were initially processed using internal Perl programs.

Bioinformatic analysis

The gene modeling annotation files and the reference genome (Oar v3.1) were accessed from Ensembl (https://asia.ensembl.org/index.html). Bowtie2 (LangmeadandS and Salzberg, 2012) was used to build the index of the reference genome, and HISAT2 (Pertea et al., 2016) was used to align paired-end clean reads to the reference genome. The mapped readings of each muscle sample were put together using StringTie (Pertea et al., 2016) with a reference-based method. Then, to identify lncRNAs, we assessed the assembled transcripts according to five criteria. The following transcripts were eliminated: 1) those with exon numbers <2, 2) those with lengths ≤200 bp, 3) those with known non-lncRNA annotations, 4) those with fragments per kilobase of exon per million fragments mapped (FPKM) < 0.5, and 5) Phylogenetic codon substitution frequency (PhyloCSF) (Sun et al., 2013). PFAM-scan (Kong et al., 2007), coding-noncoding-index (CNCI) (Finn et al., 2008), and coding potential calculator (CPC) (Kellis and Jungreis, 2011) were used to separate mRNAs from lncRNAs. All of the aforementioned tools predict that novel protein-coding transcriptome candidates are transcripts with coding potential, while novel lncRNAs are transcripts without coding potential. To determine the expression level of mRNA in each sample, FPKM (fragments per million fragments per exome) of coding genes was calculated using Cuffdiff (Trapnell et al., 2010). Differential expression analysis was performed using DESeq (Likun et al., 2010), with a default significant difference of two for fold change and a q value of 0.05.

Bowtie aligned small RNA tags to reference sequences to examine the expression and distribution of small RNA within the genome (LangmeadandS and Salzberg, 2012). Bedtools (https://bedtools.readthedocs.io/) was used to search for known miRNAs by matching them to entries in miRBase (http://www.mirbase.org/). MiRDeep2 (Friedlnder et al., 2012) was used to analyze the remaining reads to predict novel miRNAs. The following criteria were used to evaluate the expression levels of miRNAs and circRNA using TPM (Transcripts Per Kilobase of exon model per Million mapped reads) (Zhou et al., 2010). DESeq (Likun et al., 2010) was used to evaluate the differential expression, with a fold change threshold of 2 and a p value of 0.05.

CircRNAs were screened using find_circ (Memczak et al., 2013) and CIRI2 (Yuan and Zhao, 2015), and known and new circRNAs in each sample were calculated (read count) and normalized with TPM (Zhou et al., 2010). Based on the distance from the corresponding circRNA along the genome sequence, the closest protein-coding genes for circRNAs were found to mark circular RNAs. Using circBase ID, which refers to circBase annotation, all known circRNAs were given names, and novel circRNAs were given names following the rank number listed. Exonic circRNA, intronic circRNA, antisense circRNA, intergenic circRNA, 3′UTR circRNA, 5′UTR circRNA, and ncRNA circRNA were the seven categories used to categorize all circRNAs (Memczak et al., 2013). The expression levels of circRNAs were quantified by the number of junction-spanning reads, and an absolute p value < 0.05 and fold change ≥2 were considered to be significantly differentially expressed. The coding potential of circRNAs was predicted using the presence or absence of an internal ribosome entry site (IRES) independent of the 5′ cap structure on the circRNA. IRES finder (Zhao et al., 2018) was used to predict whether circRNA sequences have potential elements of IRES.

Gene ontology and kyoto encyclopedia of genes and genomes enrichment analysis

DEGs, DE circRNA parental genes, DE lncRNAs, and DE miRNA target genes were subjected to functional and pathway enrichment analysis by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GOseq was based on the Wallenius noncentral hypergeometric distribution for GO enrichment analysis (Young et al., 2010). KOBAS v3.0 (Chen et al., 2011) was used to measure the statistical significance of enrichment in KEGG pathways.

Prediction of lncRNAs, miRNAs, and CircRNAs and construction of CeRNA network

10 kb upstream and downstream of the coding genes were identified as potential regulatory targets of lncRNAs. Expression of lncRNAs used for trans-acting prediction was co-expressed with mRNA rather than correlated with mRNA location. Using Pearson’s correlation coefficients (r > 0.90 or r < −0.90) as a classifier, the expressed correlation between non-coding genes and genes was computed. Using miRanda (Enright et al., 2003), PITA (Marc et al., 2004), and RNAhybrid (Ventsislav et al., 2005), the prediction of miRNA target genes was carried out. The following threshold settings for PITA were established: maximum target length is 50,000; energy cut-off is 10, and 0.05 is the cut-off for the p-value. The following parameters were entered for RNAhybrid: utr: 3utr.fa; mir: mature. fa. To ensure the accuracy of our findings, we set the Smith-Waterman hybridization alignment match score higher than 140 and the minimum free energy of the duplex structure lower than −10. We created lncRNA/circRNA-miRNA‒mRNA networks with lncRNA/circRNA acting as decoys, miRNA acting as the core, and mRNA acting as the target. Using Cytoscape (Michaele et al., 2011), the lncRNA, circRNA, miRNA, and mRNA interactions were created and visualized.

Validation of differentially expressed transcripts

RT-qPCR validated the expression pattern of randomly selected differentially expressed transcripts from the longissimus dorsi and biceps femoris. FastKing gDNA Dispelling RT SuperMix and SuperReal PreMix Plus (Tiangen, China) mRNA detection assays were used for miRNA tests. The miRcute Plus miRNA First-strand cDNA Kit and miRcute Plus miRNA qPCR Kit (Tiangen, China) were used for microRNA tests. RT primers and a bespoke qRT‒PCR quantitative kit (GenePharma, China) were used for circRNA detection. The RNase R- and control groups (without RNase R) were produced and translated into cDNA, and each circRNA was amplified using the primers indicated to visualize the back-splicing junction of the circRNA (Myriam and Gorospe, 2018; Heng et al., 2019). The back-splicing sites were then detected by sequencing the products. GAPDH and U6 small nuclear RNA genes were chosen as endogenous control genes (all primers are shown in Supplementary Table S12). Three biological replicates and triple reactions for each sample were used for all qPCR validations. Following amplification, the products were validated by agarose gel electrophoresis and Sanger sequencing, and the 2^−ΔΔCt method was used to determine the relative transcript abundance.

Cell culture and vector construction

Primary sheep myoblasts were isolated and cultured from embryo sheep leg muscle. Myoblasts and HEK293T cells (ATCC, United States) were cultured in high-glucose DMEM supplemented with fetal bovine serum (Hyclone, United States) and double antibiotics (1% penicillin and streptomycin). To induce myoblasts differentiation, cells were switched to a differentiation medium (DMEM; 2% horse serum) in nearly 90% confluence. Small interfering RNA (siRNA) against sheep circ_0017336 were designed and synthesized by Genepharma (Suzhou, China), and a nonspecific duplex was used as the control. Sheep myoblasts were transfected with 100 nM siRNA using Lipofectamine 2000. The si-circ_0017336 sequences are 5′-UUU​GGA​GAU​AGC​AGG​GCU​GTT-3′, 3′CAG​CCC​UGC​UAU​CUC​CAA​ATT-5′. The fragment of the FBXL5 3′ UTR, including the binding site of miR-23a, was amplified and inserted into the psicheck2 vector (Promega, United States) at the 3′ end of the Renilla gene using restriction enzymes XhoI and NotI and T4 DNA ligase (TaKaRa, China). Psicheck2-FBXL5-M and psicheck2-FBXL5-W constructs were verified by sequencing.

Dual-luciferase reporter analysis

When the cell confluence reached about 80%, the miR-23a mimics and psicheck2-FBXL5-W or psicheck2-FBXL5-M were co-transfected into HEK293T cells using Lipofectamine 2000. After incubation for 24 h, Dual-luciferase activity was measured using an automatic microplate reader (Molecular Devices, United States), and the Renilla luciferase activity was normalized against firefly luciferase activity. The miR-23a mimics sequences are 5′-AUC​ACA​UUG​CCA​GGG​AUU​UCC​A-3′, 3′-GAA​AUC​CCU​GGC​AAU​GUG​AUU​U-5′.

Analytical statistics

The means and standard error of the means are used to express data (SEM). Levene’s test was used to check for homogeneity of variances, and then Student’s t-test was used to determine significance. At p < 0.05, differences were deemed statistically significant.

Results

The proportion of slow muscle fibers in the Biceps femoris and Longissimus dorsi in sheep

According to the results of immunofluorescence staining, the proportions of slow-twitch fibers in the longissimus dorsi and biceps femoris of Tan sheep were 11.4% and 37.7%, respectively, and its proportions were 5.6% and 26.8% in the longissimus dorsi and biceps femoris of Dorper sheep. Tan sheep had considerably more slow-twitch muscle fibers than Dorper sheep did (p < 0.01) in their longissimus dorsi and biceps femoris (Figures 1A,B).

FIGURE 1

FIGURE 1. The difference of the slow-twitch fiber ratio of longissimus dorsi and Biceps femoris muscle between Tan sheep and Dorper sheep. (A) Frozen sections were immunofluorescently stained with MyH 17 antibody. “DAPI” is a staining solution suitable for the nuclear staining of common cells and tissues. “Merge” is to fuse the stained nucleus with the target protein. (B) Statistical chart of immunofluorescence staining results.

Blast analysis of transcriptome sequencing

The whole-transcriptome profiling of longissimus dorsi and biceps femoris tissues from Dorper sheep and Tan sheep were obtained to evaluate the genes involved in the formation of muscle fibers. The 12 samples used for the RNA sequencing libraries yielded an average of 95.74 million clean reads, and 83.55% to 90.07% of these reads were specifically aligned to the reference genome Oar v3.1. At least 85% of the readings were equal to or greater than Q30 (Supplementary Table S1-1). A total of 22,823 mRNAs, 224 known lncRNAs, and 22,373 unique lncRNAs were found in the muscle samples after additional filtering (Figure 2A) and removal of possible coding transcripts that were identified using CNCI, CPC, and PFAM (Figure 2B). The gene structure, expression, and sequence conservation of lncRNAs and protein-coding genes were compared. As lncRNA genomic characterizations were compared to those of mRNAs, it was discovered that the length range of their transcripts was similar; more lncRNAs had 2–4 exons than mRNAs, and they also had shorter open reading frames (ORFs) and lower FPKM values (Figures 3A,B).

FIGURE 2

FIGURE 2. Identification of long non-coding RNAs (IncRNA) in muscle transcriptome. (A) Workflow for IncRNA identification. (B) Candidate IncRNAs were identified by using three applications: CNCI, coding-non-coding-index; CPC, coding potential calculator; PFAM-scan which detect and remove putative protein-coding transcripts.

FIGURE 3

FIGURE 3. Genomic architecture and expression level of IncRNAs and miRNAs. (A) Comparison of mRNA transcript length, exon number and open reading frame (ORF) length. (B) Expression level of IncRNAs and mRNAs, calculated as log10(FPKM+1). (C) Length distribution of all identified small RNAS.

With at least two independent reads crossing the back splice junction, a total of 3651 circRNAs were found. Exons, introns, and intergenic regions all serve as various circRNA source locations. Exons were more likely to produce circRNAs in this study; 3141 circRNAs were produced from exons of 1810 genes, whereas 197 circRNAs were produced from intronic regions of 175 genes, and our data indicated that the remaining circRNAs were intergenic. CircRNAs ranged in length from 192 nucleotides to 93,278 nucleotides and were mainly located on chromosome 1 (Supplementary Table S2). Interestingly, 4611 internal ribosome entry site (IRES) were predicted for circRNAs and that the highest score may reach 0.997 (Supplementary Table S12).

Additionally, for the small RNA-Seq libraries, rigorous filtering yielded an average of 16.02 million clean reads, with approximately 95% of these clean reads matched to the sheep reference sequence (Supplementary Table S1-2). Most clean reads ranged from 20 to 24 nt (Figure 3C). There were 238 known miRNAs in each sample, according to statistical statistics (Supplementary Table S2). A total of 141 mature miRNAs were annotated, and 97 new mature miRNAs were found (Supplementary Table S2).

Differential expression of genes and noncoding RNAs (LncRNAs, miRNAs, and circRNAs) between Tan sheep and Dorper sheep

In biceps femoris tissues, twelve DEGs were obtained between Tan sheep and Dorper sheep, where 11 were upregulated and 24 were downregulated (Supplementary Table S3; Figure 4A). Several DEGs were specifically expressed in Tan sheep, such as SFSWAP, or Dorper sheep, such as ASAH1. These genes may regulate the formation of muscle fibers.

FIGURE 4

FIGURE 4. Differentially expressed mRNAs and non-coding RNAs in Tan sheep and Dorper sheep. (A) Expression profiles of mRNA in biceps femoris of Tan sheep and Dorper sheep. (B) Expression profiles of IncRNA in biceps femoris of Tan sheep and Dorper sheep. (C) Expression profiles of miRNA in biceps femoris of Tan sheep and Dorper sheep. (D) Expression profiles of circRNA in biceps femoris of Tan sheep and Dorper sheep. In the volcano plots, green points represent downregulated RNAs, while red points represent upregulated RNAS. “D_B” is the biceps femoris group of the Dorper sheep. “T_B” is the biceps femoris group of the Tan sheep.

Similarly, 172 DELs were discovered in the biceps femoris between Tan and Dorper sheep, of which 28 were upregulated and 144 were downregulated (Supplementary Table S5; Figure 4B). Several DE lncRNAs were particularly expressed in Tan sheep, including LNC 000295, LNC 000300, and LNC 000306, or Dorper sheep, including LNC_017854, LNC_02118,6, and LNC_015846. In the biceps femoris, there was a total of 15 differential expression miRNAs, of which 5 were upregulated and 7 were downregulated between Tan sheep and Dorper sheep (Supplementary Table S8; Figure 4C). MiR-133 (Martin et al., 2015), miR-370 (Zhang et al., 2021), and miR-148 (Yin et al., 2020) have been revealed to be important in the development of muscle fibers among these differentially expressed miRNAs. Additionally, 95 DECs were discovered in the biceps femoris between Tan sheep and Dorper sheep, of which 59 were upregulated and 36 were downregulated (Supplementary Table S10; Figure 4D).

In longissimus dorsi tissues, twenty-five DEGs were found between Tan sheep and Dorper sheep, including 12 upregulated and 13 downregulated in Tan sheep, which have more than 2-fold differential expression (Supplementary Table S3). These DEGs included F-box and leucine-rich repeat protein 5 (FBXL5) and myosin light chain 3 (MYL3), which are related to the regulation of the formation of muscle fibers. Several DEGs were specifically expressed in Tan sheep, such as FBXL5, or Dorper sheep, such as SFSWAP.

In all, 214 DELs were obtained in the longissimus dorsi tissues of Tan sheep and Dorper sheep, of which Tan sheep had 101 upregulated and 113 downregulated DELs (Supplementary Table S5). Several DE lncRNAs were specifically expressed in Tan sheep, such as LNC_003872, LNC_01870,8, and LNC_021269, or Dorper sheep, such as LNC_008982, LNC_00889,6, and LNC_001693. These lncRNAs may control how muscle fibers are formed. Three known miRNAs were expressed overall at lower levels in Tan sheep than in Dorper sheep, whereas one known miRNA was expressed at higher levels (Supplementary Table S8). One of these miRNAs with variable expression, miR-23a, has been linked to the development of muscle fibers (Chen et al., 2021). In the longissimus dorsi tissues of Tan sheep and Dorper sheep, a total of 91 DECs were also discovered, of which 45 were upregulated and 46 were downregulated (Supplementary Table S10).

Functional analysis of differentially expressed transcripts

GO and KEGG pathway analysis were carried out to assess the potential roles of differentially expressed genes, lncRNAs, miRNAs, and circRNAs. Three functional categories were carried out of the GO pathways for DEGs (biological process, cellular component, and molecular function). The biological processes involved in the development of muscle fibers, such as calcium ion transport (GO: 0006816), oxidation-reduction process (GO: 0055114), muscular system process (GO: 0003012), and regulation of myoblast differentiation (GO: 0045661), were also significantly enriched. The CAMP signaling pathway and metabolic pathways were significantly enriched according to KEGG enrichment analysis (Supplementary Table S4; Figures 5A,B). The development of muscle fibers also includes several important metabolic and myoblast differentiation-related genes, including ATP6V0A1 and MYL3.

FIGURE 5

FIGURE 5. GO and KEGG enrichment analysis in biceps femoris of Tan sheep and Dorper sheep. (A) GO enrichment annotation of mRNAs. (B) KEGG annotation of mRNAs. “D_B”is the biceps femoris group of the Dorper sheep. “T_B” is the biceps femoris group of the Tan sheep.

Two separate algorithms—cis (genomic location) and trans (expression correlation)—were used to predict the target genes of the identified lncRNAs to elucidate their potential roles in the formation of muscle fibers, 2,154 cis-acting lncRNAs target a total of 720 genes among them, while 446 trans-acting lncRNAs have a total of 163 target genes. The GO pathway of cis-lncRNAs was highly enriched in negative regulation of glycolysis (GO: 0006110), muscle fiber development (GO: 0048741), and ATPase activity (GO: 0042623), and skeletal muscle contraction (GO: 0003009) (Supplementary Table S6). KEGG enrichment was significantly enriched in the metabolic pathway muscle atrophy; GO annotation of trans-lncRNAs was significantly enriched in positive regulation of skeletal muscle tissue development (GO: 0048643) and muscle cell differentiation (GO: 0042692); KEGG enrichment was significantly enriched in glycolysis, TCA cycle, and Wnt signaling pathway (Supplementary Table S7). GO annotation and KEGG enrichment analysis of the targets showed miRNAs were mainly enriched in metabolic pathways (Supplementary Table S9). The parental genes of DE circRNAs were largely enriched in myofibrils (GO: 0030016) and glycogen metabolic pathways (GO: 0005977), according to GO analysis. DE circRNAs were primarily enriched in the Ca²⁺ signaling pathways, FoxO signaling pathways, and AMPK signaling pathways, according to KEGG pathway analysis (Supplementary Table S11).

Validation of differentially expressed transcript expression

To validate the RNA-seq results, MYL3, ASAH1, miR-409-5p, circ_0017336, circ_0007039, LNC_0017854, LNC_003716, and LNC_000300 were selected, and their expression patterns in biceps femoris tissue of Tan sheep and Dorper sheep were examined using qPCR. FBXL5, miR-23a, miR-26a, and circ_0017430 were also selected, and their expression patterns in the longissimus dorsi of Tan sheep and Dorper sheep were examined using qPCR (Figure 6A). RT–qPCR verified that the expression patterns of randomly above differentially expressed transcripts were highly consistent with those obtained by RNA sequencing. Notably, RNase R digestion assays further demonstrated that circ_0017336 and circ_0017430 had a circular structure (Figure 6B), and the back splice junctions of circRNAs were confirmed before validation (Figure 6C).

FIGURE 6

FIGURE 6. (A) The expression levels of mRNA, miRNA, IncRNA and circRNA were verified by qPCR. (B) The results obtained from use of Rnase R indicate that the expression levels of circ_0017336 and circ_0017430 are resistant to Rnase R. “Rnase R (+)” indicates RNA treated with RNase R: “Rnase R (−)” indicates untreated RNA. (C) The backsplice junctions of circ_0017336 and circ_0017430 were verified with Sanger sequencing. The red arrows indicate head-to-tail back-splicing sites of circRNAs.

Construction of CeRNA regulatory network

We found lncRNA and circRNA gene pairs with miRNA binding sites based on the ceRNA hypothesis, and we constructed lncRNA‒miRNA-gene pairs and circRNA-miRNA-genes with lncRNA and circRNA acting as decoys, miRNA as the core, and mRNA as the target. The ceRNA network generated above was then visualized with Cytoscape software, using various shapes to represent various types of RNAs. The specific circRNA-miRNA-gene networks and lncRNA‒miRNA-gene networks closely related to the formation of muscle fibers in Tan sheep are shown in Figures 7A,B and a summary table listing all possible functioning ceRNA networks was added to Supplementary Table S14. Furthermore, from the ceRNA network, we observed some ceRNA subnetworks, which showed that circ_0017336 and its target FBXL5, ACACB and EXOC6 “talked” to each other through the same miR-23a response elements, whereas LNC_014172 and LNC_003716 “talked” to their targets through miR-23a response elements, respectively. Therefore, we speculate that these subnetworks may play a key role in the formation of muscle fiber.

FIGURE 7

FIGURE 7. (A) CircRNA-miRNA-mRNA interaction network specifically expressed in Tan sheep. (B) LncRNA-miRNA-mRNA interaction network specifically expressed in Tan sheep. The ovals, triangles, parallelograms and V types represent mRNAs, miRNAs, IncRNAs and circRNAs, respectively.

In addition, we also constructed ceRNA networks in the longissimus dorsi and biceps femoris tissues. In the longissimus dorsi, the network shows possible interactions among 2 DE lncRNAs, 3 DE miRNAs, 4 DE circRNAs, and 4 DE mRNAs (Figure 8A; Supplementary Table S14). For example, circ_0017336 were identified as ceRNA of miR-23a, which targeted MYL3, FBXL5, and C1QTNF7; LNC_003716 and LNC_014172 regulated ENSOART00000001318 by competing miRNA response elements of miR-23a. The network depicts potential connections among 4 DE lncRNAs, 7 DE miRNAs, 16 DE circRNAs, and 13 DE mRNAs in the biceps femoris (Figure 8B; Supplementary Table S14). For instance, ASAH1 and circ_0000481 potentially create a ceRNA network through miR-409-5p as a bridge. These findings imply that genes and noncoding genes that exhibit differential expression in Tan sheep and Dorper sheep might work together to control the formation of muscle fibers via an interaction network, and their regulatory patterns are also different in different tissue sites.

FIGURE 8

FIGURE 8. ceRNA network construction in differnt muscle tissues. (A) LncRNA/circRNA-miRNA-mRNA interaction network was constructed based on different expression levels in the longissimus dorsi of Tan sheep and Dorper sheep. (B) LncRNA/circRNA-miRNA-mRNA interaction network was constructed based on different expression levels in the biceps femoris of Tan sheep and Dorper sheep. The ovals, triangles, parallelograms and V types represent mRNAs, miRNAs, IncRNAs and circRNAs respectively.

miR-23a directly bound FBXL5

To investigate whether circ_0017336 acts as a competing endogenous RNA in sheep myoblasts, the interfering fragment si-circ_0017336 was transfected into sheep myoblasts to significantly inhibit the expression level of circ_0017336. The changes of miR-23a and FBXL5 expression level were detected after 3 days of differentiation of myoblasts treated with si-circ_0017336, and it was found that the interference of circ_0017336 significantly increased the expression of miR-23a in skeletal muscle cells, but suppressed the expression of FBXL5 (Figure 9A). To further validate the targeting relationship between miR-23a and FBXL5, the prediction software RNAhybrid revealed that FBXL5 had putative miR-23a binding sites (Figure 9B). To verify whether miR-23a directly targets FBXL5, miR-23a mimics were co-transfected with psi-CHECK2 dual-luciferase reporters containing the 3′UTR of FBXL5 into 293T cells for luciferase activity analysis, and we found that overexpression of miR-23a reduced the luciferase activity (Figure 9C). These results revealed that FBXL5 could indeed bind miR-23a directly.

FIGURE 9

FIGURE 9. (A) si-circ_0017336 promotes the expression of miR-23a and inhibits the exprsion of FBXLS in sheep myoblasts. (B). The binding sites of miR-23a and FBXLS were predicted by RNAhybrid software. (C) Dual-luciferase reporter gene experiments verified that miR-23a directly targets FBXL5. indicates p < 0.05, ***indicates p < 0.001, NC is the control group.

Discussion

The formation of muscle fibers is a polygenic trait in sheep, mainly determined by genetic factors, and is an important determinant of meat quality characteristics. Previous studies have revealed a closely association of the type of muscle fiber with meat quality traits due to the effects on the postmortem metabolic rate in the conversion of muscle to meat. The association between more slow-twitch fibers and better meat quality has been observed in a variety of species, including sheep (Valin et al., 1982), pigs (Bowker et al., 2005; Zhou et al., 2021), chickens (Ouyang et al., 2017), and calves (Maltin et al., 1998). Two mechanisms have been suggested for the muscle fiber pattern in Tan sheep. According to the research mentioned above, one theory proposes that the difference in the amount of slow-twitch fiber, with Tan sheep having more and Dorper sheep having less, is what causes the phenotypic. An alternative explanation claims that the difference is the result of an epigenetic regulatory process, possibly involving noncoding RNAs.

Given that some essential genes linked to specific phenotypes or significant biological processes have been identified using high-throughput RNA sequencing, it may be possible to understand the mechanisms underlying the formation of muscle fiber patterns. Numerous muscle tissue transcription profiles in animals have been studied using RNA sequencing, but similar studies have not been conducted in Tan sheep. To determine the critical elements involved in the formation of muscle fibers, we evaluated, for the first time, the differences in the expression profiles of mRNAs, lncRNAs, and circRNAs in the longissimus dorsi and biceps femoris tissues of two sheep breeds. Additionally, we found that a more complex ceRNA regulatory network is formed in the longissimus dorsi than in the biceps femoris. As studies have shown that slow-twitch fibers have higher mitochondrial content and oxidative phosphorylation capacity, and the ratio of slow-twitch fibers in the longissimus dorsi muscle is higher than that in the biceps femoris muscle, which makes it possible that the longissimus dorsi muscle has more glycolysis, oxidative phosphorylation, ad some metabolism-related genes to form a more complex ceRNA regulatory network.

Different metabolic and contractile characteristics exist between slow-twitch and fast-twitch fibers. Fast-twitch fibers have low oxidative metabolism and low mitochondrial content, while slow-twitch fibers have a high oxidative capability and large mitochondrial content. FBXL5, a member of the F-box protein family, regulates iron homeostasis in cells and systems by promoting the degradation of iron regulatory protein 2, and is an essential sensor for bioavailable iron (Ruiz and Bruick, 2014; Moroishi et al., 2014). The lack of iron in the body is closely related to the quality and function of skeletal muscle (Scherbakov et al., 2019), and we think that the formation of muscle fibers may be affected by the regulation of iron content by FBXL5. The acetyl-CoA carboxylase beta (ACACB) and ATPase H+ transporting V0 subunit a1 (ATP6V0A1) genes were significantly upregulated in Tan sheep, both of which were enriched in the oxidative phosphorylation pathway and consistent with the fact that Chinese native sheep have a better flavor than sheep from other regions. The main function of EFHB is to bind calcium ions (Rosado, 2021), triggering a Ca²⁺-dependent transport pathway (Wu et al., 2002; Chin, 2005), causing the formation of skeletal muscle fibers. Additionally, MYL3 showed differential expression in the muscle tissues of Tan sheep and Dorper sheep, which was consistent with earlier research that also deemed these genes to be key candidates for muscle growth and development (Ouyang et al., 2017). According to GO analysis, MYL3 is mostly involved in muscle contraction. Additionally, the specific roles of SFSWAP and two novel transcripts (ENSOART00000001243 and ENSOART00000001318) are unknown and are differentially expressed in both tissues. While ENSOART00000001243 and ENSOART00000001318 are downregulated in both longissimus dorsi and biceps femoris tissues, SFSWAP was interestingly increased in longissimus dorsi tissues but downregulated in biceps femoris tissues. The varied functions of these genes during the formation of these two tissues may be related to these consistent or inconsistent trend changes in expression trends, which are worthy of future analysis.

Noncoding RNAs act as epigenetic regulators of the expression of protein-coding genes in eukaryotes, which can control the expression at both transcriptional and post-transcriptional stages (Kornienko et al., 2013), and the large-scale ceRNA networks composed of noncoding RNAs and mRNA are important for regulation in various physiological and pathological processes (Salmena et al., 2011; Ala et al., 2013). Numerous miRNAs are important regulators of the growth of muscle fibers. For instance, miR-23a affects the expression of the MYHC gene and the development of muscle cells (Wang et al., 2012; Chen et al., 2021). The constructed ceRNA network results also demonstrated that DE circRNAs competed with miR-23a for binding and may function as ceRNAs to regulate the levels of ACACB, FBXL5, and EXOC6, respectively. This confirms the significance of miR-23a not only in its functions but also as a bridge in the ceRNA mechanism for the formation of muscle fibers. We hypothesized that several subnetworks within the network, such as circ_0017336-miR-23a-FBXL5, may be extremely important in controlling the formation of muscle fibers. RT-qPCR results showed that there was a ceRNA relationship between circ_0017336 and miR-23a as well as FBXL5 during sheep myogenic differentiation, and the dual luciferase reporter systems further verified that miR-23a directly targets FBXL5 for regulation. However, their roles and relationships still require confirmation. Our knowledge of the formation of muscle fiber phenotypes in sheep is expected to progress with the identification of new functions for lncRNAs and circRNAs.

Taken together, our data provide evidence for interactions that have high functional specificity in the formation of muscle fibers and that are consistent with the ceRNA hypothesis. This study provides new insights into the complex molecular mechanisms underpinning sheep meat traits variation but understanding exactly how the different type of muscle fiber are formed in the Tan sheep and Dorper sheep will require additional experiments.

Conclusion

The current study paints a thorough picture of the differences between Tan sheep and Dorper sheep in the whole-transcriptome profiles of the longissimus dorsi and biceps femoris tissues. Tan sheep had much more slow-twitch fibers than Dorper sheep had in their longissimus dorsi and biceps femoris. Among two sheep breeds, 60 DEGs were discovered, including ACACB, ATP6V0A1, ASAH1, EFHB, MYL3, C1QTNF7, SFSWAP, and FBXL5. Act according to our research, numerous lncRNAs, miRNAs, and circRNAs play a crucial role in crucial biological processes related to the development of muscle fiber pathways, such as the FoxO signaling pathway, the AMPK signaling pathway, and the Ca²⁺ signaling pathway. We also created the first lncRNA/circRNA-miRNA-gene interaction network based on the differentially expressed transcripts from muscle tissues, where some networks, such as circ_0017336-miR-23a-FBXL5, may be crucial in the control of muscle fiber creation. Our findings have established future investigations on the molecular mechanisms underlying muscle fiber development by identifying possible regulators and molecular regulatory networks that may be connected to the production of muscle fibers in sheep.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The animal study was reviewed and approved by animal welfare committee of the State Key Laboratory for Agro-biotechnology of China Agricultural University (approval number XK257).

Author contributions

XK, YuL, XL, YaL and DF carried out the tests; RC, SC, YW, and ZL carried out the data analysis; ML, FL, and MF conceptualized the study and supervised the experiments. RC, YW, and ZL also created the figures and tables. All authors reviewed and approved the final version of the manuscript.

Funding

The National Key Research and Development Program of China (2021YFD1200900) and the Program of Agricultural Breeding in the Ningxia Hui Autonomous Region provided financial assistance for this work (No. nxnyyz20150101).

Conflict of interest

The reviewer QZ declared a shared affiliation with the author(s) YuL to the handling editor at time of review.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.991606/full#supplementary-material

Abbreviations

ceRNA, competing endogenous RNA; DE, endogenoustially expressed genes; DEMs, Differentially expressed miRNAs; DELs, Differentially expressed lncRNAs; DECs, Differentially expressed circRNAs; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; qPCR, quantitative Polymerase Chain Reaction; ACACB, acetyl-CoA carboxylase beta; ATP6V0A1, ATPase H+ transporting V0 subunit a1; ASAH1, N-acylsphingosine amidohydrolase 1; EFHB, EF-hand domain family member B; MYL3, myosin light chain 3; FBXL5, F-box and leucine rileucine-richotein 5; SFSWAP, splicing factor SWAP; C1QTNF7, C1q and TNF relateTNF-rerelatedong non-coding RNAs; mRNAs, messenger RNAs; miRNAs, microRNAs; circRNAs, circular RNA; MyH7, myosin heavy chain 7; MyHC, myosin heavy chain, cardiac muscle complex.

References

Ala, U., Karreth, F. A., Bosia, C., Pagnani, A., Taulli, R., Leopold, V., et al. (2013). Integrated transcriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments. Proc. Natl. Acad. Sci. U. S. A. 110 (18), 7154–7159. doi:10.1073/pnas.1222509110

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowker, B. C., Swartz, D. R., Grant, A. L., and Gerrard, D. (2005). Myosin heavy chain isoform composition influences the susceptibility of actin-activated S1 ATPase and myofibrillar ATPase to pH inactivation. Meat Sci. 71, 342–350. doi:10.1016/j.meatsci.2005.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, B., Ma, M., Zhou, Z., Kong, S., Zhang, J., Zhang, X., et al. (2022). circPTPN4regulates myogenesis via the miR-499-3p/NAMPTaxis. J. Anim. Sci. Biotechnol. 13 (1), 2. doi:10.1186/s40104-021-00664-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). Kobas 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, 316–322. doi:10.1093/nar/gkr483

CrossRef Full Text | Google Scholar

Chen, M., Wei, X., Song, M., Jiang, R., Huang, K., Deng, Y., et al. (2021). Circular RNA circMYBPC1 promotes skeletal muscle differentiation by targeting MyHC. Mol. Ther. Nucleic Acids 24, 352–368. doi:10.1016/j.omtn.2021.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chin, E. R. (2005). Role of Ca2+/calmodulin-dependent kinases in skeletal muscle plasticity. J. Appl. Physiol. 99 (2), 414–423. doi:10.1152/japplphysiol.00015.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Doerr, C. R. AshmoreandL. (1971). Comparative aspects of muscle fiber types in different species. Exp. Neurol. 31 (3), 408–418. doi:10.1016/0014-4886(71)90243-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5 (11), R1. doi:10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Finn, R. D., Tate, J., Mistry, J., Coggill, P. C., Sammut, S. J., Hotz, H. R., et al. (2008). The Pfam protein families database. Nucleic Acids Res. 32 (1), D281–D288. doi:10.1093/nar/gkm960

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedlnder, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40 (1), 37–52. doi:10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Heng, L., Yin, H., Jing, Y., Yan, X., Chen, W., Jiang, C., et al. (2019). Profiles analysis reveals circular RNAs involving zebrafish physiological development. J. Cell. Physiol. 234 (9), 15922–15933. doi:10.1002/jcp.28250

CrossRef Full Text | Google Scholar

Hong, L., Gu, T., He, Y., Zhou, C., Hu, Q., Wang, X., et al. (2019). Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development. Front. Cell Dev. Biol. 7 (289), 289. doi:10.3389/fcell.2019.00289

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, J., Lv, W., Xia, P., Xu, Z. Y., Zheng, A. D., Wang, X. J., et al. (2018). Long noncoding RNA SYISL regulates myogenesis by interacting with polycomb repressive complex 2. Proc. Natl. Acad. Sci. U. S. A. 115 (42), E9802–E9811. doi:10.1073/pnas.1801471115

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz, J. C., and Bruick, K. (2014). F-box and leucine-rich repeat protein 5 (FBXL5): Sensing intracellular iron and oxygen. J. Inorg. Biochem. 133, 73–77. doi:10.1016/j.jinorgbio.2014.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, X., Liu, G., Liu, Y., Xu, Q., Zhang, M., and Fang, M. (2013). Transcriptome profile at different physiological stages reveals potential mode for curly fleece in Chinese tan sheep. PLoS One 8 (8), e71763. doi:10.1371/journal.pone.0071763

PubMed Abstract | CrossRef Full Text | Google Scholar

Kellis, M., and Jungreis, I. (2011). PhyloCSF: A comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics 27 (13), i275–i282. doi:10.1093/bioinformatics/btr209

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, L., Zhang, Y., Ye, Z., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi:10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornienko, A. E., Guenzl, P. M., Barlow, D. P., and Pauler, F. M. (2013). Gene regulation by the act of long non-coding RNA transcription. BMC Biol. 11 (2), 59–350. doi:10.1186/1741-7007-11-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotb, A., Amaresh, P., Supriyo, D., Grammatikakis, I., Kim, J., Ding, J., et al. (2015). Circular RNAs in monkey muscle: Age-dependent changes. Aging 7 (11), 903–910. doi:10.18632/aging.100834

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristensen, L. S., Andersen, M. S., Stagsted, L., Ebbesen, K. K., Hansen, T. B., and Kjems, J. (2019). The biogenesis, biology and characterization of circular RNAs. Nat. Rev. Genet. 20 (7), 675–691. doi:10.1038/s41576-019-0158-7

PubMed Abstract | CrossRef Full Text | Google Scholar

LangmeadandS, B., and Salzberg, L. (2012). - Fast gapped-read alignment with Bowtie 2. Nat. Methods 9 (4), 357–359. doi:10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Larzul, C., Roy, P. L., Gogué, J., Talmant, A., Jacquet, B., Lefaucheur, L., et al. (1999). Selection for reduced muscle glycolytic potential in Large White pigs. II. Correlated responses in meat quality and muscle compositional traits. Genet. Sel. Evol. 31 (1), 61. doi:10.1186/1297-9686-31-1-61

CrossRef Full Text | Google Scholar

Likun, W., Zhixing, F., Xi, W., Wang, X., and Zhang, X. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinforma. Oxf. Engl. 26 (1), 136–138. doi:10.1093/bioinformatics/btp612

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Zhang, J., Xu, Q., Kang, X., Wang, K., Wu, K., et al. (2018). Integrated miRNA-mRNA analysis reveals regulatory pathways underlying the curly fleece trait in Chinese tan sheep. BMC Genomics 19 (1), 360. doi:10.1186/s12864-018-4736-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltin, C. A., SinclairandP, K. D., Warriss, D., Grant, C. M., Porter, A. D., Delday, M. I., et al. (1998). The effects of age at slaughter, genotype and finishing system on the biochemical properties, muscle fibre type characteristics and eating quality of bull beef from suckled calves. Anim. Sci. 66 (2), 341–348. doi:10.1017/s1357729800009462

CrossRef Full Text | Google Scholar

Marc, R., Peter, S., Matthias, H., and Giegerich, R. (2004). Fast and effective prediction of microRNA/target duplexes. RNA (New York, NY) 10 (10), 1507–1517. doi:10.1261/rna.5248604

CrossRef Full Text | Google Scholar

Martin, H., Jan, N., and Julie, B.-V. (2015). Muscle-specific microRNAs in skeletal muscle development. Dev. Biol. 410 (1), 1–13. doi:10.1016/j.ydbio.2015.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495 (7441), 333–338. doi:10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaele, S., Keiichiro, O., Johannes, R., Wang, P. L., and Ideker, T. (2011). Cytoscape 2.8: New features for data integration and network visualization. Bioinforma. Oxf. Engl. 27 (3), 431–432. doi:10.1093/bioinformatics/btq675

CrossRef Full Text | Google Scholar

Moroishi, T., Yamauchi, T., Nishiyama, M., and Nakayama, K. I. (2014). HERC2 targets the iron regulator FBXL5 for degradation and modulates iron metabolism. J. Biol. Chem. 289, 16430–16441. doi:10.1074/jbc.M113.541490

PubMed Abstract | CrossRef Full Text | Google Scholar

Myriam, P. A. C., and Gorospe, M. (2018). Detection and analysis of circular RNAs by RT-PCR. Bio. Protoc. 8 (6), e2775. doi:10.21769/BioProtoc.2775

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, H., Wang, Z., Chen, X., Yu, J., Li, Z., and Nie, Q. (2017). Proteomic analysis of chicken skeletal muscle during embryonic development. Front. Physiol. 8, 281. doi:10.3389/fphys.2017.00281

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11 (9), 1650–1667. doi:10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosado, J. A. (2021). SARAF and EFHB modulate store-operated Ca 2+ entry and are required for cell proliferation, migration and viability in breast cancer cells. Cancers 13, 4160. doi:10.3390/cancers13164160

PubMed Abstract | CrossRef Full Text | Google Scholar

RyuandB, Y. C., and Kim, C. (2005). The relationship between muscle fiber characteristics, postmortem metabolic rate, and meat quality of pig longissimus dorsi muscle. Meat Sci. 71 (2), 351–357. doi:10.1016/j.meatsci.2005.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: The rosetta stone of a hidden RNA language? Cell 146 (3), 353–358. doi:10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Scherbakov, N., Sandek, A., Valentova, M., Mayer, A., von Haehling, S., Jankowska, E., et al. (2022). Iron deficiency and reduced muscle strength in patients with acute and chronic ischemic stroke. J. Clin. Med. 11 (3), 595. doi:10.3390/jcm11030595

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, L., Li, H., Xin, X., Duan, Y., Hua, X. Q., and Jin, Y. (2013). Muscle fiber types, characteristics and meat quality. Adv. Mat. Res. 2203 (634-638), 1263–1267. doi:10.4028/www.scientific.net/amr.634-638.1263

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41 (17), e166. doi:10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28 (5), 511–515. doi:10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Valin, C., Touraille, C., Vigneron, P., and Ashmore, C. R. (1982). Prediction of lamb meat quality traits based on muscle biopsy fibre typing. Meat Sci. 6 (4), 257–263. doi:10.1016/0309-1740(82)90036-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ventsislav, R., Vesselin, B., and Nikiforov, M. I. (2005). MicroInspector: A web tool for detection of miRNA binding sites in an RNA sequence. Nucleic acids Res. 33, W696–W700. doi:10.1093/nar/gki364

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Chen, X., Zheng, Y., Li, F., Lu, Z., Chen, C., et al. (2012). MiR-23a inhibits myogenic differentiation through down regulation of fast myosin heavy chain isoforms. Exp. Cell Res. 318 (18), 2324–2334. doi:10.1016/j.yexcr.2012.06.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Kanatous, S. B., Thurmond, F. A., Gallardo, T., Isotani, E., Bassel-Duby, R., et al. (2002). Regulation of mitochondrial biogenesis in skeletal muscle by CaMK. Science 296 (5566), 349–352. doi:10.1126/science.1071163

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, H., He, H., Cao, X., Shen, X., Han, S., Cui, C., et al. (2020). MiR-148a-3p regulates skeletal muscle satellite cell differentiation and apoptosis via the PI3K/AKT signaling pathway by targeting Meox2. Front. Genet. 11 (4), 512. doi:10.3389/fgene.2020.00512

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol. 11 (2), R14. doi:10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, G., and Zhao, J. (2015). Ciri: An efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16 (1), 4. doi:10.1186/s13059-014-0571-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P., Du, J., Guo, X., Wu, S., He, J., Li, X., et al. (2021). LncMyoD promotes skeletal myogenesis and regulates skeletal muscle fiber-type composition by sponging miR-370-3p. Genes 12 (4), 589. doi:10.3390/genes12040589

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J., Wu, J., Xu, T., Yang, Q., He, J., and Song, X. (2018). IRESfinder:Identifying RNA internal ribosome entry site in eukaryotic cell using framed k-mer features. J. Genet. Genomics 45 (7), 403–406. doi:10.1016/j.jgg.2018.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Chen, J., Li, Z., Li, X., Hu, X., Huang, Y., et al. (2010). Integrated profiling of MicroRNAs and mRNAs: MicroRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One 5 (12), e15224. doi:10.1371/journal.pone.0015224

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Liu, Y., Zhang, L., Kong, X., and Li, F. (2021). Serine-to-glycine ratios in low-protein diets regulate intramuscular fat by affecting lipid metabolism and myofiber type transition in the skeletal muscle of growing-finishing pigs. Anim. Nutr. 7 (7), 384–392. doi:10.1016/j.aninu.2020.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar