Comprehensive Analysis of miRNA and mRNA Expression Profiles during Muscle Development of the Longest Dorsal Muscle in Gannan Yaks and Jeryaks

1. Introduction

The Jeryak is the result of crossbreeding between Gannan yak and Jersey cattle. Numerous studies have demonstrated that from Jersey cattle for crossbreeding with Gannan yak, the hybrid offspring of Jeryaks exhibit significantly increased body size and weight, robust physique, strong adaptability, rapid growth and development [1], which provides support for herders’ living production and pastoral economy [2]. In view of this, exploring the molecular regulatory mechanisms of skeletal muscle development may provide a direction to improve the production performance of Jeryaks.

The economic value of mammals is evaluated based on the meat production and its quality, which is influenced by the growth of skeletal muscles [3]. The growth of skeletal muscle is an intricate and accurate process, including the creation of fetal muscle fibers prior to delivery and the enlargement of muscle fibers after delivery [4]. Furthermore, the regulation of numerous transcription factors is imperative for the growth of skeletal muscles [5]. Muscle derived regulatory factors (MRFs) and muscle cell enhancer 2 (MEF2) are widely recognized as crucial elements in the growth and development of skeletal muscle, as evidenced by previous studies [6,7,8]. And PAX3 and PAX7 act as regulatory factors that control the development of skeletal muscle by interacting with MRFs [9]. Muscle growth inhibitor (MSTN) has been identified as a specific negative regulatory factor involved in skeletal muscle development, which interacts with certain signaling pathways to inhibit the expression of MRFs or PAX [10,11]. Furthermore, the Wnt, MAPK, and PI3K-Akt signaling pathways play a pivotal role in myogenesis, regeneration and differentiation across multiple species including mice, chicken and humans [12,13,14,15,16,17]. In recent years, the advancements in high-throughput sequencing technology have revealed that skeletal muscle growth is influenced not only by protein-coding RNAs but also by non-coding RNAs (ncRNAs), specifically microRNAs (miRNAs). These ncRNAs play a significant role in various biological processes and are regulated through multiple mechanisms.

MiRNAs, a class of small RNAs that are highly conserved across evolutionary species, exert post-transcriptional regulation on target genes by either degrading or inhibiting the translation of specific messenger RNAs (mRNAs) [18,19]. It has been demonstrated to play a role in vital biological processes, including cell proliferation and differentiation [20,21,22].Currently, a plethora of research studies have indicated that the abundant expression of muscle-specific miRNAs, such as miR-1, miR-206, miR-133, and miR-499, within skeletal muscles significantly influences their development [23,24,25]. Chen et al. [26] reported an observed notable increase in the levels of miR-1 and miR-206 during satellite cells differentiation in skeletal muscle, while a decrease was observed during the process of skeletal muscle regeneration. Further investigation unveiled that miR-1/206 exerted inhibitory effects on the expression of Pax7, thereby hindering the proliferation and differentiation of satellite cells. Nakasa et al.[27] showed that the activation of MyoD, MyoG and Pax7 occurred when miR-1, miR-133, and miR-206 were injected together at the site of muscle injury in rats and this activation subsequently facilitated muscle regeneration and suppressed muscle fibrosis. Meanwhile, MiR-499 has been demonstrated its significance in regulating the composition of muscle fibers [28]. Additionally, miR-22-3p impedes the proliferation of skeletal muscle cells and promotes cellular differentiation by virtue of its interaction with IGFBP3 [29]. The identified interactions between oar-mir-655-3p and oar-mir-381-5p with ACSM3 and ABAT play a crucial role in the process of muscle organogenesis and myoblast migration in sheep [30]. In recent years, transcriptome sequencing has been extensively explored to investigate the correlation between miRNAs and mRNAs in skeletal muscle development. However, the regulatory mechanisms governing miRNA-mRNA interaction in relation to Jeryaks muscle growth and development have yet to be investigated.

Therefore, in this study, we generated six small RNA libraries by collecting samples from the longissimus dorsi muscle of adult Gannan yaks and Jeryaks. By employing high-throughput sequencing and conducting comprehensive bioinformatics analysis, we elucidated the expression patterns of miRNA and mRNA, ultimately constructing intricate miRNA-mRNA regulatory networks between miRNA and mRNAs that are relevant to muscle growth. These findings help us to deepen our understanding of the hybrid dominance of Jeryaks and provide a theoretical basis for further research on the regulatory mechanisms of miRNAs associated with Jeryak muscle growth and development.

2. Materials and Methods

2.1. Sample Collection

Three adult Jeryaks and three Gannan yaks (P1, P2, P3, M1, M2, and M3) were selected from the Qinaihe Female Cattle-yak Breeding Base, Hezuo City, Gannan Tibetan Autonomous Prefecture, Gansu Province, and divided into Groups P and M, with the same environmental conditions and rearing conditions. The selected experimental animals were slaughtered according to the slaughterhouse’s procedure, and the longissimus dorsi muscle were collected on the slaughtering line. All samples were gathered in 1.5 mL tubes for centrifugation, promptly transferred to liquid nitrogen, and later preserved at -80 °C for subsequent sequencing.

2.2. Construction and Sequencing of Small Ribonucleic acid Library

Following the utilization of the Trizol kit (Invitrogen, Carlsbad, CA, USA) to extract the complete RNA from the samples, the polyacrylamide gel electrophoresis (PAGE) technique was employed to enrich the bands within the 18-30nt range. After connecting the 3’ and 5’ connectors, the small RNA linked to both connectors underwent reverse transcription and PCR amplification. In the end, a group of competent bands of about 140 bp was obtained and purified after gelation with PAGE, then dissolved in an EB solution. Following the construction of the library, the samples underwent testing for eligibility utilizing the Agilent2100 and ABI StepOnePlus Real-Time PCR systems provided by Life Technologies. In the end, the sequencing was carried out utilizing the Illumina Hiseq2000 machine (Illumina, San Diego, CA, USA).

2.3. Differentially Expressed miRNA and mRNA Analysis

We will process the information obtained from sequencing. To begin with, eliminate reads that have multiple bases of poor quality (q-value ≤ 20) or unknown nucleotides (N); secondly, the part of the reads that do not contain the 3’ aptamer and the part that contain the 5’ aptamer should be removed; subsequently, the reads that contain the 3’ and the 5’ linkers without small RNA fragments in the middle should be removed; then, reads containing ployA should be removed; and lastly, reads with sizes less than 18 nt should be removed. To confirm and eliminate non-coding RNAs, such as rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs, we compared all clean labels with small RNAs in the databases through the GenBank (v209.0) and Rfam (v11.0) databases. All cleaning markers were aligned with the yak reference genome (BosGru_v3.0). Meanwhile, labels that were mapped to exons, introns, or repeating sequences were eliminated. Then，novel miRNAs were predicted by miRDeep2 software [31]. And calculated its miRNA expression level, normalized to per million transcripts (TPM). We analyzed (DE) miRNAs and DEMs in Jeryak and Gannan yak using DESeq2 [32] software. In the comparison, we identified miRNAs with |fold change| > 2.0 and P < 0.05 as significant (DE) miRNAs. Meanwhile, we defined mRNAs with |fold change| > 2.0 and FDR < 0.05 as significant DEMs.

2.4. Target Gene Prediction of Differentially Expressed miRNAs and KEGG and GO Analysis

We utilized both Miranda (v3.3a) [33] and TargetScan (v7.0) [34] to conduct target gene prediction for miRNAs that were already known, identified, and newly discovered. We considered the intersections of the prediction results as target genes of miRNAs. Next, the GO (http //www.geneontology.org/) and KEGG (http //www.genome.jp/kegg/) databases were utilized to ascertain significant biological functions and pathways. And with a screening threshold of P < 0.05, the screened terms and pathways were defined as biological processes that were significantly enriched in the target genes of differentially expressed miRNAs.

2.5. miRNA–mRNA Interaction Network Construction

We screened miRNA-mRNA targeting relationship pairs by co-analysis of miRNA and mRNA and screened muscle growth-related mRNAs based on literature search. miRNA expression correlation with target genes was then evaluated using Pearson correlation coefficient (PCC). If the expression of a target gene is inversely related to the expression of the corresponding miRNA, we consider it as a potential target gene for differentially expressed (DE)miRNAs. Using Cytoscape 3.1.0 (http //www.cytoscape.org/), the mapping of interaction networks between miRNA and mRNA was finally accomplished.

2.6. Real-Time qPCR Validation of Differentially Expressed mRNA and miRNAs

In order to confirm the accuracy of the sequencing findings, we conducted qRT-PCR validation on a random sample of 15 genes that exhibited differential expression, consisting of 8 mRNAs and 7miRNAs. Initially, we employed Primer 5.0 to design primers for mRNAs. And the forward primers for miRNAs were designed by miRNA Design V1.01. The primers for DEMs and (DE) miRNAs are listed in Table 1 and Table 2, respectively. Afterwards, total RNA was extracted from the samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). MiRNA reverse transcription kits (Transgen Biology, Beijing, China) were used to reverse transcribe and fluorescently quantify the extracted total RNA. The reactions were performed with U6 as an internal reference for miRNA and GAPDH as an internal reference for mRNA. Furthermore, we analyzed miRNA expression levels in different tissues to confirm that miRNAs are involved in bovine skeletal muscle development. Briefly, we isolated total RNA from heart, spleen, lung, liver, kidney, longestissimus dorsi muscle and subcutaneous fat using TRIzol. Subsequently, we converted it to cDNA by reverse transcription and finally evaluated miRNA expression using qRT-PCR with U6 as internal reference. Each reaction was repeated three times. Graphs were generated for all trials using GraphPad Prism 9, and data were analyzed using the 2-^ΔΔCT method.

3. Results

3.1. Summary of Sequencing Small RNA

To detect (DE)miRNAs during muscle development in Gannan yaks and Jeryaks, we created and sequenced six small RNA libraries (M1, M2, M3, P1, P2, and P3) using the Illumina HiSeq™ 2000 platform. A total of 17,408,914, 14,015,050, 16,597,591, 14,318,748, 9,195,586, and 16,740,804 clean readings were produced in M1, M2, M3, P1, P2, and P3, correspondingly. Following the elimination of adapters, impurities, and reads of poor quality, a total of 16255077, 13362108, 15172291, 13782038, 8756687, and 16180410 clean tags were acquired and utilized for subsequent examination (Table 3). After comparing all clean tags with GenBank and Rfam, it was discovered that a total of 74809937 (89.61%) known miRNAs and 19881 (0.02%) novel miRNAs were obtained from Gannan yaks and Jeryaks once small RNAs matching rRNA, tRNA, snRNA, snoRNA, and scRNA were eliminated (Figure 1A, Table S4). In the meantime, the yak reference genome (BosGru_v3.0) had a comparable match to 75.23% of the cleaning markers. The size distribution of small RNAs in the longissimus dorsi muscle of two breeds is comparable, with the majority being concentrated around 22 nt (Figure 1B). To further determine the quality of sequencing data, pairwise Pearson correlations were calculated between any pair of small RNA seq samples, and the results showed similarity between the same variety while differences between the two varieties (Figure 1C). Principal component analysis (PCA) revealed similarities between samples and differences between two cattle breeds (Figure 1D). The presence of these findings suggests that the sequencing data is of excellent quality.

3.2. Differential Expression Analysis of mRNAs

Group M and group P had a total of 1,819 genes that were expressed differentially (FDR < 0.05). Among them, compared with group M, 738 genes were up-regulated in group P, while the remaining 1081 genes were down-regulated in group P (Figure 2A). The results of GO enrichment analysis for the DEMs are shown in (Figure 2B). 1819 DEMs are annotated in 987 GO terms (P < 0.05) (Table S5). Among them, 77 GO terms are significantly enriched in cell components, including terms such as cells, cell parts, and smooth muscle contraction fibers. In the category of molecular functions, 91 molecular functional terms have been enriched, mainly including insulin like growth factor I binding, protein activity, protein binding, and other terms. A total of 819 GO terms showed significant enrichment in biological processes, primarily related to metabolic processes, development of skeletal muscle tissue, migration of muscle cells, differentiation of skeletal muscle cells, proliferation and differentiation of skeletal muscle cells, contraction of smooth muscles, and proliferation of muscle cells. By utilizing the KEGG database, the analysis of molecular pathways for genes with differential expression revealed that a noteworthy total of 20 signaling pathways exhibited significant enrichment (P < 0.05) (Figure 2C, Table S6). The genes of differential expression primarily engaged in the MAPK signaling pathway, Regulation of actin cytoskeleton, and Apelin signaling pathway.

3.3. Differential Expression Analysis of miRNAs

We compared miRNAs from groups M and P using DESeq2(V1.20.0) software and identified 230 differentially expressed miRNAs, consisting of 172 already characterized miRNAs and 58 newly discovered ones. In comparison to Gannan yaks, the longissimus dorsi muscle of Jeryaks have 135 miRNAs were down-regulated and 95 miRNAs were up-regulated (Figure 3A). And they are clustered in the two branches of M and P respectively (Figure 3B).

3.4. Target Gene Prediction of Differentially Expressed miRNAs and KEGG and GO Analysis

As is well known, miRNA controls biological processes by suppressing the translation of target genes [35]. Hence, to gain a deeper comprehension of the possible functions of miRNAs in the growth of skeletal muscle, Targetscan and miRanda were employed for the anticipation of its target genes. Contained a sum of 17285 target genes that were predicted by the 230 (DE)miRNAs (Table S7). The analysis of GO enrichment indicates that these target genes are associated with 3217 functional categories that have been significantly enriched (P < 0.05) (Table S8). Among them, 270 GO terms were enriched in the cellular component classification, including cell, cell part, Wnt signalosome, PRC1 complex, striated muscle thin filament and other terms. Within the category of molecular function, there were 575 enriched terms, mainly including catalytic activity, binding to growth factors, protein interaction, activity as a transcription regulator, binding to transcription factors, binding to Ras GTPase, and other related terms. There is a significant enrichment of 2372 GO terms in biological processes, primarily involving the promotion of cellular process, development of multicellular organisms, control of gene expression, developmental processes, Wnt-mediated cell-cell signaling, development of animal organs, and other related terms (Figure 4A). The analysis of KEGG indicated that the target genes of (DE)miRNA showed a notable enrichment in 118 signaling pathways (P < 0.05) (Table S9).The enriched pathways mainly involved Ras signaling pathway, MAPK signaling pathway, Wnt signaling pathway, FoxO signaling pathway and Fatty acid metabolism(Figure 4B).

3.5. Validation of differentially expressed miRNAs and mRNAs by qRT-PCR

To establish the validity of the RNA-Seq outcomes, we evaluated seven differentially expressed miRNAs, including 4 up-regulated miRNAs (miR-450-x, miR-1271-z, miR-136-x, miR-142-y) and 3 down-regulated miRNAs (miR-339-x, miR-98-y, miR-204-x). As well as 8 differentially expressed mRNA, including 4 up-regulated mRNA (IGFBP5, FOSB, MYBPH, ANKRD1) and 4 down-regulated mRNA (CFL2, DDX23, DDX47, UBE2G2). The results showed that the miRNA and mRNA expression patterns of Gannan yak and Jeryak were similar to the sequencing results (Figure 5A). This result further supports the repeatability and reliability of our sequencing data. And then We investigated the expression of miR-2478-z by qRT-PCR, which targets multiple genes involved in muscle growth. Our study found that miR-2478-z expression was relatively high in the longest dorsal muscle tissue (Figure 5B). Thus, these also results confirm the reliability of the reciprocal network and suggest that (DE)miRNAs and DEMs play a regulatory role in muscle growth and development.

3.6. Building the Network for miRNA-mRNA Interactions

These 230 (DE) miRNAs were analyzed in association with previous DEMs based on the negative regulatory effects of miRNAs on mRNAs. To further explore the potential role of miRNAs in Jeryak muscle growth. There are 5916 pairs of miRNA-mRNA relationships have been identified (Table S10). To gain a more extensive comprehension of the interplay between miRNAs and target genes, we chose 20 differentially expressed mRNAs (DEMs) associated with muscle development and 14 (DE)miRNAs that target them based on a thorough examination of literature. The miRNA-mRNA interaction network was constructed using this information (Figure 6). It should be emphasized that some miRNAs, like miR-2478-z, miR-339-x, along with some newly discovered miRNAs (novel-m0036-3p and novel-m0037-3p), have the ability to target numerous genes associated with muscle development. The findings suggest that these miRNAs could have a significant impact on the regulation of muscle development.

4. Discussion

Muscle growth is a multifaceted economic characteristic that not only impacts the quality of meat from livestock and poultry, but also plays a pivotal role in enhancing overall meat production performance[36,37]. Compared to conventional bovines, yaks exhibit a slower growth rate and have relatively modest production performance[38]. Therefore, the quest for strategies to optimize yak productivity continues to be a pressing inquiry in need of resolution. Jeryak, a hybrid breed resulting from crossbreed of Gannan yak and Jersey cattle, have exhibited enhanced production capabilities in comparison to Gannan yak[39]. Hence, investigating the molecular regulatory mechanism underlying this phenomenon will contribute to enhancing the economic advantages associated with yaks.

As widely acknowledged, the process of muscle growth is intricately governed by pivotal genes and their associated signaling pathways[40,41]. A total of 1819 DEMs were identified (FDR < 0.05) in this study. Despite the limited sample size of this study, a plethora of genes with different molecular functions were identified, thereby providing a valuable reference for investigating the signaling pathways associated with the longest dorsal muscle in Jeryak. In this study, a number of GO terms and KEGG pathways have been identified as closely associated with the process of muscle growth Examples include skeletal muscle cell differentiation, skeletal muscle tissue development and, MAPK and the signaling pathway. Notably, the background genes associated with these pathways have also been found to be involved in muscle growth and development of muscles, including MEF2C, SIX1, CFL2, and FGF6. It has been shown that the MEF2C gene plays an active role in skeletal muscle differentiation and regeneration [42,43]. Yang et al. [44]. demonstrated that knockdown of the MEF2C gene inhibited myogenic differentiation of myoblasts, and conversely, overexpression of MEF2C promoted their differentiation. In addition, the Six1 gene, a potent regulator of skeletal muscle development [45], activates myogenesis and enhances Myf3 expression by specifically binding to the MEF5 site in the promoter of the gene, thereby promoting skeletal muscle development [46,47]. Sun et al. [48]. reported a significant upregulation of CFL2 expression in skeletal muscle, which plays a crucial in regulating the activity of transcription factors involved in muscle formation within C2C12 cells and this, in turn, has an impact on the proliferation and differentiation of myofibroblasts. Another study revealed a correlation between the inactivation of FGF6 and impaired muscle damage repair, as well as reduced expression of MyoD and MyoG in mutant mice and it was suggested that this could be attributed to the inhibition of satellite cell proliferation and differentiation caused by the inactivation of FGF6 [49]. Although these studies have contributed to prediction of key gene the function and accuracy, further investigation is required to explore the regulation mechanisms of other functional differentially expressed mRNAs (DEMs) associated with muscle growth and development in the longest dorsal muscle.

MiRNAs, crucial non-coding regulators in post-transcriptional regulation, exert their influence on a plethora of biological processes including cell proliferation, apoptosis, and the development of tumors [50]. A total of 230 miRNAs exhibited differential expression in the present study. Among these, a total of 135 (DE)miRNAs displayed noteworthy down-regulation in the longestissimus dorsi muscle of Jeryak, whereas 95 (DE)miRNAs were exhibited significant up-regulated. Performing GO and KEGG pathway enrichment analysis on differentially expressed miRNA target genes to investigate their potential functional roles. The GO enrichment results demonstrate that (DE)miRNA primarily participates in the regulation of Wnt signaling and developmental processes. Additionally, KEGG pathway analysis revealed significantly enrichment of some (DE)miRNA in key pathways including Ras signaling, MAPK signaling, Wnt signaling, and FoxO signaling. As a prominent pathway for miRNA target genes in Jeryak, the Ras signaling cascade has been documented to exert inhibitory effects on skeletal muscle myogenesis [51]. The MAPK signaling pathway has the capacity to regulate biological processes via multiple mechanisms of activating or inhibiting associated factors [52]. It has been found that activation of the p38/MAPK signaling pathway can contribute to the thickening of muscle fiber cross-section and the increase of muscle length by regulating the protein content in adult muscle fibers after birth, thereby increasing total muscle mass [53,54]. A previous study demonstrated that the Wnt signaling exerts direct control over the expression of myogenic regulatory factors (MRFs) during l embryonic development in animals, thereby influencing muscle production [55]. Additionally, the FoXO signaling pathway is also involved in muscle growth and development, and it leads to skeletal muscle atrophy primarily through protein degradation [56]. Sandri et al. [57]. demonstrated that FoXO induces skeletal muscle atrophies by upregulating the ubiquitin ligase Atrogin-1.

To further scrutinize the candidate miRNAs governing muscle growth and development in Jeryak, we successfully constructed a miRNA-mRNA interaction network based on the results of our comprehensive analysis of mRNAs and miRNAs. Significantly, the two newly identified miRNAs exhibited differential expression in both species. The function of their target genes may mediate the effects of these miRNAs on muscle growth and development. Within the intricate network of interactions, novel-m0036-3p and novel-m0037-3p were found to specifically target both SIX2, SIX1, NDRG4, STX4, and CTCF genes. SIX1 has been reported as a constituent of the vertebrate SIX gene family. During vertebrate skeletal myogenesis and development, the Six1 transcription factor exerts a transcriptional regulatory effect on the myogenic determinant gene family, thereby indirectly modulating downstream muscle development-related genes through its regulation of individual members of MRFs [58]. Recent investigation reveals that excessive expression of the transcription factor SIX2 stimulates satellite cell growth in bovine skeletal muscles [59]. Zhu et al. [60]. reported that NDRG2 treatment of C12C4 myoblasts activated the Akt/CREB signaling pathway and a significantly upregulated MyoD and MyoG gene expression, thereby promoting myoblast differentiation. STX4, also known as Syntaxin 4, exhibits significant expression levels in the skeletal muscle tissue. Yoo et al. [60] found that STX4 enhanced the proliferation and differentiation of myoblasts by interacting with CDO. CTCF is critical for early embryonic development [61]. CTCF has been found to be a factor involved in the regulation of myogenesis, which promotes muscle differentiation mainly by coordinating with MRFs [62]. Interestingly, we also identified SIX2, SIX1, NDRG4, STX4, and CTCF as down-regulated DEMs. As a result, we speculated that these miRNAs might influence the growth and development of muscles by inhibiting the expression of DEMs associated with muscles.

In addition, related studies have shown that not only one miRNA can act on multiple mRNAs, but also multiple miRNAs can act on the same Mrna [63]. Notably, our interaction network diagram showed that miR-339-x, miR-339-z, and miR-10926-z acted concurrently with MICAL2, and miR-421-y and miR-450-x acted concurrently with MEF2A. Therefore, we hypothesize that these miRNAs may bind specifically to target mRNAs and have an effect on Jeryak skeletal muscle growth and development. Future studies should involve validation of possible miRNA-mRNA targeting relationship pairs in subsequent studies to explore their mechanisms of action in skeletal muscle growth and development.

5. Conclusions

In this study, a total of 230 (DE)miRNAs and 1819 DEMs were identified by RNA-Seq. The results of both GO and KEGG analyses have shown that (DE)miRNAs are principally engaged in signaling pathways related to muscle growth and development. Subsequently, 14 key miRNAs were screened by constructing miRNA-mRNA interaction networks, including miR-2478-z, miR-339-x, novel-m0036-3p and novel-m0037-3p. These miRNAs may regulate muscle growth and development in Jeryaks by targeting their target genes.