Comparative Transcriptome Analysis of Slow-Twitch and Fast-Twitch Muscles in Dezhou Donkeys

The skeletal muscle fiber profile is closely related to livestock meat quality. However, the molecular mechanisms determining muscle fiber types in donkeys are not completely understood. In this study, we selected the psoas major muscle (PM; mainly composed of oxidative-type muscle fibers) and biceps femoris muscle (BF; mainly composed of glycolytic-type muscle fibers) and systematically compared their mRNA and microRNA transcriptomes via RNA-seq. We identified a total of 2881 differentially expressed genes (DEGs) and 21 known differentially expressed miRNAs (DEmiRs). Furthermore, functional enrichment analysis showed that the DEGs were mainly involved in energy metabolism and actin cytoskeleton regulation. The glycolysis/gluconeogenesis pathway (including up-regulated genes such as PKM, LDHA, PGK1 and ALDOA) was more highly enriched in BF, whereas the oxidative phosphorylation pathway and cardiac muscle contraction (including down-regulated genes such as LDHB, ATP2A2, myosin-7 (MYH7), TNNC1, TPM3 and TNNI1) was more enriched in PM. Additionally, we identified several candidate miRNA–mRNA pairs that might regulate muscle fiber types using the integrated miRNA–mRNA analysis. Combined with the results of protein–protein interaction (PPI) analysis, some interesting DEGs (including ACTN3, TNNT3, TPM2, TNNC2, PKM, TNNC1 and TNNI1) might be potential candidate target genes involved in the miRNA-mediated regulation of the myofibril composition. This study is the first to indicate that DEmiRs, especially eca-miR-193a-5p and eca-miR-370, and potential candidate target genes that are mainly involved in actin binding (e.g., ACTN3, TNNT3 and TNNC1) and the glycolysis/gluconeogenesis pathways (e.g., PKM) might coregulate the myofibril composition in donkeys. This study may provide useful information for improving meat quality traits in Dezhou donkeys.


Introduction
Skeletal muscle consists of different types of muscle fiber, which may affect the biochemical characteristics of the muscle [1]. The major muscle fibers can be roughly divided into the slow-twitch type (type I) and fast-twitch type (type II) according to the differences in their metabolic and contractile characteristics. With regard to metabolism, slow-twitch fibers, which are also known as the oxidative muscle type, show a higher abundance of mitochondria and myoglobin and higher oxidative metabolism than fast-twitch fibers (glycolytic type) [2]. It is well established that muscle fiber composition is an important factor affecting meat quality, including meat color, tenderness and flavor [3,4]. In addition to meat quality in animals, muscle fiber composition is also strongly related to muscle health in animals [5]. Therefore, the control of muscle fiber composition is of special concern not only for meat quality regulation, but also for molecular diagnostics and therapeutics. It was reported that gene expression in eukaryotes is specific to each tissue [6]. Further, the amount of gene products that are made in the same tissue, as well as in other tissues

Animal Selection and Sample Collection
The trial was carried out at the Dezhou donkey original breeding farm authorized by Shandong Province (Dezhou city, Shandong, China) and all procedures used in this study were approved by the Institutional Animal Care Committee at Liaocheng University (Permit No. DFG21010103-1). Six healthy male donkeys with similar body weight (225 ± 6 kg; mean ± SEM) and aged 2.5 years were selected for this study. After a 12 h feed withdrawal, six donkeys were slaughtered for tissue collection, as was done in our previous publication [21]. After death, the neck muscle (NM), psoas major muscle (PM), longissimus doris muscle (LD) and biceps femoris muscle (BF) were sampled and washed with sterilized saline, and immediately placed in liquid nitrogen for histochemical analysis and RNA extraction. Then, all frozen samples stored in dry ice were transported to the laboratory and stored at −80 • C for further analysis.

Histochemical Analysis
Histochemical analyses were performed with ATPase staining using the method of Brooke and Kaiser [22] with slight modifications. Transverse serial sections of 10 µm in thickness were cut from entire blocks (1.0 × 1.0 × 1.5 cm) with a cryostat microtome (HM525, Microm GmbH, Walldorf, Germany) at −20 • C. Sections were subsequently used for histochemical analysis of myosin adenosine triphosphatase (mATPase) following alkaline (pH 10.70) preincubation. An image analysis system (Image-Pro ® plus 5.1, Media Cybernetics Inc., Rockville, MD, USA) was used to examine the stained sections. The muscle fibers were classified as fiber types I, IIA and IIB according to the nomenclature of Brooke and Kaiser. Approximately 500 fibers per sample were counted to analyze the muscle fiber characteristics with Image-Pro Plus software (Image-Pro ® plus 5.1, Media Cybernetics Inc., Rockville, MD, USA), including the fiber number percentage (FNP) and fiber area percentages (FAP). FNP shows the ratio of the counted fiber number of each fiber type to the total counted fiber number. FAP was the ratio of a total cross-sectional area of each fiber type to total fiber area measured.

RNA Extraction, Sequencing and Bioinformatics Analysis
Total RNA was extracted from PM and BF muscle samples of the donkeys using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. RNA quality was verified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and the NanoDrop2000 (Thermo Scientific, Wilmington, NC, USA), and then high-quality RNA samples were used for library preparation and sequencing. The paired-end RNA-seq sequencing libraries were sequenced with the Illumina NovaSeq 6000 (2 × 150 bp read length) using 5 µg of total RNA. The RNA-seq data were analyzed with the free online platform Majorbio Cloud Platform (http://www.majorbio.com). Briefly, quantification of gene expression as transcript per million (TPM) values was carried out using the RSEM algorithm (Version 1.3.1) [23]. The differential expression genes (DEGs) were identified and screened by DESeq2 based on fold change > 1.5 and p-value < 0.05 [24].

Functional Enrichment Analysis
All DEGs and target genes of DEmiRs were used for functional enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, with the free online platform Majorbio Cloud Platform (http://www.majorbio.com). The results with p-adjust values less than 0.05 were considered significantly enriched.

Real-time Quantitative PCR
To validate the transcriptome results, six mRNAs and two miRNAs were randomly selected for qRT-PCR. The relative expression analysis of genes was assessed by qRT-PCR as described in our previous study [9]. Trizol reagent was used to extract the total RNA of muscle samples, and reverse transcription was carried out using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China). GAPDH gene was used as a housekeeping gene in this procedure to normalize the gene expression data. Real-time qPCR was performed using ChamQ SYBR Color qPCR Master Mix (2X) (Vazyme, Nanjing, China) in a 7300 Real-Time PCR System (Applied Biosystems, USA). For miRNA analysis, 1 µg of RNA was used for reverse transcription with Maxima Reverse Transcriptase (Thermo Fisher Scientific, Vilnius, Lithuania). The U6 snRNA was used as the internal reference. Real-time qPCR was performed using 2X SG Fast qPCR Master Mix (Bimake, Houston, TX, USA) in a QuantStudio 3&5 Real-Time PCR System (ABI/Thermo Fisher, Wilmington, NC, USA). The primers used in this study are listed in Table S1. The 2 -∆∆Ct method was used to quantitate mRNA and miRNA expression level.

Protein-Protein Interaction (PPI) Network Analysis
To further identify hub genes, a protein-protein interaction (PPI) network analysis was performed in the Search Tool for the Retrieval of Interacting Genes (STRING) database. In the network, each edge was obtained based on the edge weight score to quantify the interaction confidence, and the PPI network was constructed and visualized using Cytoscape software (v3.9.0, Cytoscape Consortium, San Diego, CA, USA).

Statistical Analyses
SPSS statistical software (Version 22, SPSS IBM, New York, NY, USA) was used for statistical analysis. Data on the muscle fiber characteristics of NM, LD, PM and BF muscle was assessed by one-way ANOVA followed by Duncan multiple comparison, and Student's t test was used to analyze the differences in FAP and RT-PCR data between the PM and BF. The variability of results was expressed as the mean ± standard error. Means were considered significantly different at p < 0.05.

Data Availability
The dataset(s) supporting the conclusions of this article are available in the National Center for Biotechnology Information Sequence Read Archive (http://www.ncbi.nlm. nih.gov/sra) under BioProject accession numbers PRJNA870101 (RNA-seq data) and PR-JNA870031 (small RNA-seq data).

Muscle Fiber Type Populations
Muscle fiber characteristics, including fiber diameter (FD), fiber number percentage (FNP) and fiber area percentage (FAP) are summarized in Table S2. In total, 24 samples (six sets of four muscles) were collected and analyzed, including the neck muscle, longissimus dorsi muscle from the back, psoas major muscle from the lower levels of the lumbar spine, and biceps femoris muscle from the hind limbs. As expected, a clear difference in the muscle fiber type composition was observed among the four muscle sections. The type I (slow) FAP of Dezhou donkeys was most prevalent in PM (55.99%), followed by NM (34.82%), LD (21.94%) and BF (9.22%) (Table S2 and Figure 1). Additionally, the diameter of the PM was thinner than that of the BF. Based on the above results, PM and BF were selected for the further study of muscle fiber types. tions. The type I (slow) FAP of Dezhou donkeys was most prevalent in PM (55.99%), followed by NM (34.82%), LD (21.94%) and BF (9.22%) (Table S2 and Figure 1). Additionally, the diameter of the PM was thinner than that of the BF. Based on the above results, PM and BF were selected for the further study of muscle fiber types.

Overview of RNA and Small-RNA Sequencing
To assess the expression of mRNAs and miRNAs in the different muscle fiber types in donkeys, we collected PM and BF for the transcriptomic profiling of all mRNAs and miRNAs through the high-throughput sequencing of eight muscle samples (four each from PM and BF). In the mRNA-Seq libraries, a total of 504,369,474 raw reads were produced from the eight samples analyzed (Table S3) and deposited in the NCBI database under accession number PRJNA870101. After filtration, a total of 499,455,950 clean reads were obtained from the eight samples tested, more than 85% of which were mapped to the Dezhou donkey reference genome (Table S3). The Q20 and Q30 clean read quality scores of the eight samples were above 97% and 94%, respectively, indicating that the reliability and quality of the sequencing data met the standard analysis requirements and could be used for further analysis.
Additionally, in the small RNA-Seq libraries, an average of 11,452,819 clean reads were obtained from the eight samples. An average of more than 300 known miRNAs and 70 novel miRNAs were obtained after analyses and annotation (Table S4). Furthermore, the data supporting these results were deposited in the NCBI database under accession number PRJNA870031.

Overview of RNA and Small-RNA Sequencing
To assess the expression of mRNAs and miRNAs in the different muscle fiber types in donkeys, we collected PM and BF for the transcriptomic profiling of all mRNAs and miRNAs through the high-throughput sequencing of eight muscle samples (four each from PM and BF). In the mRNA-Seq libraries, a total of 504,369,474 raw reads were produced from the eight samples analyzed (Table S3) and deposited in the NCBI database under accession number PRJNA870101. After filtration, a total of 499,455,950 clean reads were obtained from the eight samples tested, more than 85% of which were mapped to the Dezhou donkey reference genome (Table S3). The Q20 and Q30 clean read quality scores of the eight samples were above 97% and 94%, respectively, indicating that the reliability and quality of the sequencing data met the standard analysis requirements and could be used for further analysis.
Additionally, in the small RNA-Seq libraries, an average of 11,452,819 clean reads were obtained from the eight samples. An average of more than 300 known miRNAs and 70 novel miRNAs were obtained after analyses and annotation (Table S4). Furthermore, the data supporting these results were deposited in the NCBI database under accession number PRJNA870031.

Differentially Expressed Genes (DEGs) between the Two Muscles
The PCA showed that the DEGs differed between the two muscle types of donkeys ( Figure S1). Muscles were differentiated based on the presence of slow and fast myofibers. A total of 2881 differentially expressed genes (DEGs) were identified by comparing the gene expression between PM and BF, among which 1575 genes were up-regulated and 1306 genes were down-regulated in BF relative to PM (Figure 2A and Table S6a). The results indicated that there was a significant difference between the two muscle types at the mRNA level. As expected, BF showed higher expression levels of genes related to the fast-twitch muscle fiber type (up-regulated DEGs, p < 0.05), such as TNNI2, TNNC2, MYH1, MYH2, ACTN3, MYBPC2 and MYLPF, whereas PM showed higher abundances of genes related to the slow-twitch muscle fiber type (down-regulated DEGs, p < 0.05), including MYH7, MYH7B, TNNC1, TNNI1, MYL2 and MYL3.
( Figure S1). Muscles were differentiated based on the presence of slow and fast my bers. A total of 2881 differentially expressed genes (DEGs) were identified by compar the gene expression between PM and BF, among which 1575 genes were up-regula and 1306 genes were down-regulated in BF relative to PM (Figure 2A and Table S The results indicated that there was a significant difference between the two mu types at the mRNA level. As expected, BF showed higher expression levels of genes lated to the fast-twitch muscle fiber type (up-regulated DEGs, p < 0.05), such as TNN TNNC2, MYH1, MYH2, ACTN3, MYBPC2 and MYLPF, whereas PM showed hig abundances of genes related to the slow-twitch muscle fiber type (down-regula DEGs, p < 0.05), including MYH7, MYH7B, TNNC1, TNNI1, MYL2 and MYL3. To better understand the function of these DEGs, we performed further GO a KEGG pathway enrichment analyses ( Figures 2B,C and S2). After GO enrichment ana sis, we observed that the DEGs were mainly enriched in biological pathways. Amo the identified biological processes, the enriched GO terms were mainly related to mu development and energy metabolism ( Figure S2A). The results of the KEGG enrichm analysis also showed that the DEGs were enriched in some pathways involved in mu profiles and function as well as energy metabolism, such as the glucagon signal pathway, PPAR signaling pathway, glycolysis/gluconeogenesis, AMPK and PI3K-A signaling pathways, and regulation of actin cytoskeleton ( Figure S2B). Additionally, glycolysis/gluconeogenesis pathway (including up-regulated genes such as PK To better understand the function of these DEGs, we performed further GO and KEGG pathway enrichment analyses ( Figure 2B,C and Figure S2). After GO enrichment analysis, we observed that the DEGs were mainly enriched in biological pathways. Among the identified biological processes, the enriched GO terms were mainly related to muscle development and energy metabolism ( Figure S2A). The results of the KEGG enrichment analysis also showed that the DEGs were enriched in some pathways involved in muscle profiles and function as well as energy metabolism, such as the glucagon signaling pathway, PPAR signaling pathway, glycolysis/gluconeogenesis, AMPK and PI3K-AKT signaling pathways, and regulation of actin cytoskeleton ( Figure S2B). Additionally, the glycolysis/gluconeogenesis pathway (including up-regulated genes such as PKM, BPGM, ENO2, LDHA, PGK1 and ALDOA) was more highly enriched in BF, whereas the oxidative phosphorylation pathway (including down-regulated genes such as NDUFB3, NDUFC1, SDHB, ATP5MC1, NDUFB5 and ATP5F1D) was more enriched in PM ( Figure 2B,C). In addition, some important genes involved in the calcium signaling pathway (such as CAMK1, CALM1, CALML6, ANXA2 and NFATC3 genes) were also identified in our study, and most of them were more highly expressed in BF.

Differentially Expressed miRNAs (DEmiRs) and Target Gene Prediction
To understand the expression characteristics of microRNAs in the slow-twitch PM and fast-twitch BF muscles, we constructed a small RNA library of the two muscles. After sequence alignment, 2528 known miRNAs and 621 novel miRNAs were identified (Table  S4). A total of 21 known miRNAs were identified as DE-miRs between the two muscles (Table S5). Among these miRNAs, 7 miRNAs (such as eca-miR-199b-5p, eca-miR-370 and eca-miR-758) were up-regulated, and 14 miRNAs (such as eca-miR-196b, eca-miR-196a, eca-miR-192, eca-miR-615-3p, eca-miR-499-3p, eca-miR-128 and eca-miR-193-5p) were down-regulated in BF. Furthermore, 48,194 miRNA target genes were predicted and identified based on analyses with miRanda and RNAhybrid software. Among these genes, 1531 target genes were found to be previously identified DEGs, which were considered intersection genes (Table S6b).

qRT-PCR Validation
To validate the transcriptome results, the relative expression of five mRNAs (including ACTN3, CPT1B, MYL1, MYL2 and PFKL) and two miRNAs (eca-miR-10a and eca-miR-758) was randomly selected for qRT-PCR quantification. As shown in Figure S3, the qRT-PCR expression patterns of these DEGs and two DEmiRs were in accordance with the results of the RNA-Seq analysis, indicating the reproducibility and reliability of our RNA-Seq results. Notably, MYH4, a marker of fast-twitch muscle fiber, was significantly more abundant in BF than in PM ( Figure S3A), despite the failure to identify this gene in RNA-seq analysis.

Combined Analysis of DEGs and DEmiRs
To reveal possible miRNA-mRNA interactions, we used Cytoscape mapping software (v3.9.0, Cytoscape Consortium, San Diego, CA, USA) to construct the potential miRNA-mRNA regulatory network. Based on GeneCards (www.genecards.org) and some references [6,26], more than 100 DEGs that might be closely related to muscle fiber composition were screened from the mRNA sequencing data (Table S6c). By considering these genes combined with the intersection genes, we screened 74 genes to construct a miRNA-mRNA regulatory network (Table S6d). We obtained 64 pairs of DEmiRs-DEGs, which included 9 down-regulated and 2 up-regulated miRNAs as well as 16 down-regulated and 48 upregulated DEGs, and we then drafted the network of the DEmiRs and the DEGs based on the negative correlations ( Figure 3). In the networks, eca-miR-193-5p, eca-miR-1379 and eca-miR-370 showed large numbers of target genes. Some important genes of the muscle fiber profile were found in this network, such as ACTN3, TNNC2, TNNT3, MYPLF, TNNC1, NFATC2, PKM, TPM2 and TNNI1. Most genes that were up-regulated in BF were targets of the down-regulated miRNAs eca-miR-193-5p and eca-miR-1379 as well as eca-miR-615-3p ( Figure 3B). On the other hand, TNNC1, TPM2, TNNI1 and MYH14 (also referred to as MYH7b), enriched in PM, were targets of the up-regulated miRNA eca-miR-370 ( Figure 3A). 2022, 13, x FOR PEER REVIEW 8 of 14 Additionally, GO and KEGG enrichment analyses were carried out to explore the potential functions of these DEGs in this network. As shown in Figure S4A, the enriched GO terms were mainly related to biological processes including the Wnt signaling pathway, cell differentiation and glycolytic processes, cellular components such as troponin complex, and molecular functions such as actin binding. The enriched KEGG pathways were mainly associated with energy metabolism (such as glycolysis/gluconeogenesis and the pentose phosphate pathway), the Wnt signaling pathway, cancer (mainly related to cell proliferation and differentiation), and muscle contraction ( Figure S4B). Finally, these genes in the miRNA-mRNA regulatory network were also used to build a PPI network, which included 272 edges and 50 nodes. The top 10 hub genes in network string interactions ranked by the degree method were ACTN3, TPM2,TNNT3, TNNC2, MYLPF, TNNC1, TNNI1, TNNI2, PKM and MMP2 (Figure 4). The degree scores are shown in Table S6e. Additionally, GO and KEGG enrichment analyses were carried out to explore the potential functions of these DEGs in this network. As shown in Figure S4A, the enriched GO terms were mainly related to biological processes including the Wnt signaling pathway, cell differentiation and glycolytic processes, cellular components such as troponin complex, and molecular functions such as actin binding. The enriched KEGG pathways were mainly associated with energy metabolism (such as glycolysis/gluconeogenesis and the pentose phosphate pathway), the Wnt signaling pathway, cancer (mainly related to cell proliferation and differentiation), and muscle contraction ( Figure S4B). Finally, these genes in the miRNA-mRNA regulatory network were also used to build a PPI network, which included 272 edges and 50 nodes. The top 10 hub genes in network string interactions ranked by the degree method were ACTN3, TPM2, TNNT3, TNNC2, MYLPF, TNNC1, TNNI1, TNNI2, PKM and MMP2 (Figure 4). The degree scores are shown in Table S6e.

Discussion and Conclusions
Skeletal muscle fiber is not only closely related to human muscle health but also affects meat quality in livestock production [1,4,15]. In our previous studies, we observed differences in meat quality among various muscles of donkeys (e.g., donkey gluteus vs. longissimus dorsi) [20] as well as genes related to donkey meat tenderness, as reported by Chai et al. [9]. However, the molecular mechanisms determining the muscle fiber profiles of donkeys remain unclear, especially with respect to miRNA regulation. Thus, we further explored the molecular mechanisms affecting muscle fiber types by comparing microRNA and mRNA differences through RNA sequencing of the two types of muscles. First, according to the observed muscle fiber characteristics, we selected PM and BF as representatives of type I muscle fibers (oxidative type) and type II muscle fibers (glycolytic type), respectively. Similar to previous reports, PM was characterized as the slow oxidative fiber type (type I), as found in pork [27,28], beef [29] and goat muscles [30], whereas BF showed higher respective percentages of type II fibers (fast-twitch fibers) [4,13,31]. Generally, the diameter of type I muscle fibers is thinner than that of type II muscle fibers [32], and our data were consistent with this pattern. Then, we identified a total of 2881 differentially expressed genes (DEGs) and 21 known differentially expressed miRNAs (DEmiRs) as well as several candidate miRNA-mRNA pairs that might regulate muscle fiber type by integrated miRNA-mRNA analysis.
In the present study, we also investigated the high functional divergence between PM and BF based on functional enrichment analysis. As expected, the functional analysis based on KEGG results showed that DEGs involved in energy metabolism, oxidative phosphorylation and cardiac muscle contraction were more enriched in PM, including some important genes (LDHB, ATP2A2, myosin-7 (MYH7), MYH6, TNNC1, TPM3 and TNNI1) that are generally more highly expressed in slow-twitch skeletal muscle and cardiac muscle [18,33]. However, the expression of myosin-2 (MYH2), myosin-1-like (MYH1) and MYH4, which are molecular markers of adult fast/type II myofibers [34], was significantly higher in BF, indicating that BF was more enriched with type II myofibers than PM and showed a higher glycolysis level. Concerning the functional pathways of fast-twitch muscle fibers, BF was also mainly enriched in pathways related to the glucagon signaling pathway and glycolysis/gluconeogenesis, including the PKM, PGK1, LDHA, PPP3CA, TNNC2 and ALDOA genes. For example, PKM is a key rate-limiting

Discussion and Conclusions
Skeletal muscle fiber is not only closely related to human muscle health but also affects meat quality in livestock production [1,4,15]. In our previous studies, we observed differences in meat quality among various muscles of donkeys (e.g., donkey gluteus vs. longissimus dorsi) [20] as well as genes related to donkey meat tenderness, as reported by Chai et al. [9]. However, the molecular mechanisms determining the muscle fiber profiles of donkeys remain unclear, especially with respect to miRNA regulation. Thus, we further explored the molecular mechanisms affecting muscle fiber types by comparing microRNA and mRNA differences through RNA sequencing of the two types of muscles. First, according to the observed muscle fiber characteristics, we selected PM and BF as representatives of type I muscle fibers (oxidative type) and type II muscle fibers (glycolytic type), respectively. Similar to previous reports, PM was characterized as the slow oxidative fiber type (type I), as found in pork [27,28], beef [29] and goat muscles [30], whereas BF showed higher respective percentages of type II fibers (fast-twitch fibers) [4,13,31]. Generally, the diameter of type I muscle fibers is thinner than that of type II muscle fibers [32], and our data were consistent with this pattern. Then, we identified a total of 2881 differentially expressed genes (DEGs) and 21 known differentially expressed miRNAs (DEmiRs) as well as several candidate miRNA-mRNA pairs that might regulate muscle fiber type by integrated miRNA-mRNA analysis.
In the present study, we also investigated the high functional divergence between PM and BF based on functional enrichment analysis. As expected, the functional analysis based on KEGG results showed that DEGs involved in energy metabolism, oxidative phosphorylation and cardiac muscle contraction were more enriched in PM, including some important genes (LDHB, ATP2A2, myosin-7 (MYH7), MYH6, TNNC1, TPM3 and TNNI1) that are generally more highly expressed in slow-twitch skeletal muscle and cardiac muscle [18,33]. However, the expression of myosin-2 (MYH2), myosin-1-like (MYH1) and MYH4, which are molecular markers of adult fast/type II myofibers [34], was significantly higher in BF, indicating that BF was more enriched with type II myofibers than PM and showed a higher glycolysis level. Concerning the functional pathways of fast-twitch muscle fibers, BF was also mainly enriched in pathways related to the glucagon signaling pathway and glycolysis/gluconeogenesis, including the PKM, PGK1, LDHA, PPP3CA, TNNC2 and ALDOA genes. For example, PKM is a key rate-limiting enzyme that catalyzes the last step in glycolysis [35], whereas PGK is involved in the first step in ATP generation. PKM has also been reported to participate in regulating meat quality; for instance, the expression of this gene is significantly negatively correlated with the final meat pH [28]. Another gene, LDHA, encodes skeletal-muscle-specific lactate dehydrogenase, which also plays an important role in glycolysis and affects the meat quality [36]. It is worth noting that most DEGs involved in the calcium signaling pathway, including CALM1, CALM6, CALMK1, PPP3CA, ANXA2, NFATc2 and NFATc3, showed significantly higher levels in BF (mainly composed of fast-twitch muscle fibers), which was similar to the results of previous studies showing that the CaN/NFAT signaling pathway plays a key role in regulating the muscle fiber phenotype [15,37]. It has also been shown that the expression of CALM1, a subunit of CaM, is negatively correlated with intramuscular fat in pigs [38] and influences meat quality. Furthermore, genes such as MYH1, MYH7, MYH4 and MYL2 have been identified as candidate genes affecting meat tenderness in our previous work [9]. These results indicated that there is a difference between PM and BF at the transcription level, which is consistent with the data from the histochemical analysis of the study subjects. Additionally, some interesting genes may play an important role in regulating skeletal myofiber formation and thereby affect meat quality.
MicroRNAs can regulate target gene expression by promoting mRNA degradation or repressing translation. Several studies have also identified some differentially expressed miRNAs between glycolytic and oxidative muscles [13,[15][16][17]. At the miRNA level, the present study is the first to identify a total of 21 known DEmiRs in comparisons of PM (mainly composed of slow-twitch muscle fibers) and BF (mainly composed of fast-twitch muscle fibers). Similarly, in Mongolian horses, only 11 microRNAs (including miR-499-3p, miR-499-5p and miR-206) were reported as DemiRs between the splenius (higher slow-twitch muscle fiber population) and gluteus medius (higher fast-twitch muscle fiber population) muscles [18]. Indeed, some of these DEmiRs have been reported to regulate muscle development and functions, including myofiber-type switching. For instance, miR-499 can affect myoblast differentiation and skeletal muscle fiber composition by repressing target genes such as SOX6 and FNIP [16,18,39]. In turn, eca-miR-196a and eca-miR-196b were also more highly expressed in PM, whereas their expression was most absent in BF, which may be important for controlling myofiber type. In chickens, miR-196-5p is most enriched in slow-twitch muscle fibers, and can regulate myofiber type by targeting CAML1 [15]. MiR-370-3p has been shown to down-regulate ACADSB gene expression and inhibit the formation of slow-twitch muscle fiber [40]. Additionally, some of the DEmiRs identified in our study, such as eca-miR-192, eca-miR-10a and eca-miR-128, are involved in muscle cell proliferation, differentiation and energy expenditure [41][42][43]. These results suggested that some DEmiRs identified in our study might also be involved in skeletal muscle fiber development and composition in donkeys.
To further understand the regulatory mechanism of miRNAs and their targets affecting the muscle fiber phenotype, an interaction network was constructed between the miRNAs and mRNAs involved in muscle development and muscle fiber composition. In this interaction network, some genes related to myofiber type pathways associated with the troponin complex (such as TNNI2, TNNC2 and TNNT3), actin binding (ACTN3) and glycolysis/gluconeogenesis (such as GRK6, PKM, EPHA, ALDOA and PGM1), which are negatively regulated by eca-miR-193a-5p or/and eca-miR-1379, showed higher expression in BF (higher fast muscle fiber population). According to previous reports, miRNA-193a-5p mostly functions in the context of cell proliferation and differentiation, including tumor or cancer development and 3T3-L1 preadipocyte proliferation and differentiation [44][45][46]. For instance, the overexpression of miR-193a-5p can inhibit 3T3-L1 preadipocyte proliferation and differentiation by targeting ACAA2 [46]. A study by Ju et al. showed that miR-193-3p may contribute to the regulation of oxidative myofibers in chickens by targeting the PPARGC1A gene [47]. Additionally, only one report has identified miR-1379 expressed in serum as a novel equine tissue miRNA, and the expression of this miRNA was relatively low in both muscles evaluated in our study, indicating that its regulatory role in muscle fiber composition in donkeys needs to be further confirmed. On the other hand, this study indicated that only one up-regulated miRNA showed a high degree in the interaction network. We found that some interesting genes (e.g., TPM2, TNNC1, TNNI1, MYH14 and ACADS) were predicted to be targets of eca-miR-370, and these genes were more highly expressed in slow muscle fibers [18,40,48]. For instance, Zhang et al. revealed that miR-370-3p can inhibit the formation of slow-twitch myofibers by down-regulating the ACADSB gene, which might play an important role in fiber-type transitions [40]. Furthermore, combined with the results from PPI analysis, a number of potential candidate target genes affecting myofibril composition were screened, including ACTN3, TNNT3, TPM2, TNNC2, PKM, TNNC1 and TNNI1. Considering these results together, we inferred that eca-miR-193a-5p and eca-miR-370 are very likely involved in controlling muscle fiber composition by targeting genes that are mainly involved in actin binding and the glycolysis/gluconeogenesis pathways in donkeys, although the cellular functional validation of these findings and the elucidation of the complex regulatory mechanisms involved will require further study.
In summary, we first compared the mRNA and miRNA transcript profiles of slowtwitch (PM) and fast-twitch (BF) muscles by RNA sequencing and constructed a miRNA-mRNA network. We identified 2881 DEGs and 21 DEmiRs between two muscles of Dezhou donkey. Furthermore, we identified two microRNAs (eca-miR-193a-5p and eca-miR-370) and some potential candidate target genes (ACTN3, TNNT3, TPM2, TNNC2, PKM, TNNC1 and TNNI1), which were mainly involved in actin binding and the glycolysis/gluconeogenesis pathways, and might coregulate the muscle fiber types. Therefore, these miRNAs and potential candidate target genes could be used as candidate biomarkers controlling myofibril composition, and this study may expand our understanding of the molecular mechanisms underlying meat quality traits in donkeys. and PFKL) were randomly selected for qRT-PCR validation, and GAPDH was chosen as the reference gene. Additionally, the expression of MYH4, an important marker of fast-twitch muscle fiber, was also quantified by RT-qPCR, although this gene was not identified in the mRNA sequencing data in this study. (C and D) Validation of miRNA by qRT-PCR. Two differentially expressed miRNAs (eca-miR-10a and eca-miR-758) were randomly selected for qRT-PCR validation, and U6 snRNA was used as the internal reference. PM, psoas major; BF, biceps femoris muscle. Data are expressed as the means ± SEMs. a, b Means with different letters are significantly different, p < 0.05; Figure S4: Top 20 significantly enriched GO terms (A) and KEGG pathways (B) of DEGs in the miRNA-mRNA network; Table S1: List of primers used for differentially expressed genes and differentially expressed miRNA validation; Table S2: Muscle fiber characteristics of Dezhou donkey NM, LD, PM and BF muscles; Table S3: Overview of the mRNA sequencing data; Table S4: Overview of the small RNA sequencing data; Table S5: Differentially expressed miRNAs between PM and BF; Table S6