Front. Genet.Frontiers in GeneticsFront. Genet.1664-8021Frontiers Media S.A.77255710.3389/fgene.2021.772557GeneticsOriginal ResearchComprehensive Analysis of mRNA, lncRNA, circRNA, and miRNA Expression Profiles and Their ceRNA Networks in the Longissimus Dorsi Muscle of Cattle-Yak and YakHuang et al.Difference in Expression ProfilingHuangChun1GeFei1MaXiaoming1DaiRongfeng1DingkaoRenqing2ZhaxiZhuoma3BurenchaoGetu3BaoPengjia1WuXiaoyun1GuoXian1ChuMin1YanPing1*LiangChunnian1*
Key Laboratory of Yak Breeding Engineering Gansu Province, Lanzhou Institute of Husbandry and Pharmaceutical Science, Chinese Academy of Agricultural Sciences, Lanzhou, China
Livestock Institute of Gannan Tibetan Autonomous Prefecture, Hezuo, China
Haixi Agricultural and Animal Husbandry Technology Extension Service Center, Qinghai, China
Edited by:Marcos De Donato, Instituto de Tecnología y Educación Superior de Monterrey (ITESM), Mexico
Reviewed by:Ran Di, Institute of Animal Sciences (CAAS), China
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Cattle-yak, as the hybrid offspring of cattle (Bos taurus) and yak (Bos grunniens), demonstrates obvious heterosis in production performance. Male hybrid sterility has been focused on for a long time; however, the mRNAs and non-coding RNAs related to muscle development as well as their regulatory networks remain unclear. The phenotypic data showed that the production performance (i.e., body weight, withers height, body length, and chest girth) of cattle-yak was significantly better than that of the yak, and the economic benefits of the cattle-yak were higher under the same feeding conditions. Then, we detected the expression profiles of the longissimus dorsi muscle of cattle-yak and yak to systematically reveal the molecular basis using the high-throughput sequencing technology. Here, 7,126 mRNAs, 791 lncRNAs, and 1,057 circRNAs were identified to be differentially expressed between cattle-yaks and yaks in the longissimus dorsi muscle. These mRNAs, lncRNA targeted genes, and circRNA host genes were significantly enriched in myoblast differentiation and some signaling pathways related to muscle development (such as HIF-1 signaling pathway and PI3K-Akt signaling pathway). We constructed a competing endogenous RNA (ceRNA) network and found that some non-coding RNAs differentially expressed may be involved in the regulation of muscle traits. Taken together, this study may be used as a reference tool to provide the molecular basis for studying muscle development.
Yak, a special germplasm resource mainly inhabiting the Qinghai-Tibet Plateau, has been optimized for living and a source of living for the local herdsmen. The cattle-yak constitutes the hybrid of cattle (Bos taurus) and yak (Bos grunniens) exhibiting outstanding hybrid vigor in growth rate, meat performance, plateau adaptability, etc. The meat of cattle-yak is highly enriched in protein but has lower fat than yak, meeting the requirements of a popular healthy and high-quality diet (Song et al., 2019; Wang et al., 2021b). Inhabiting the high-altitude environment, cattle-yak provides a natural green food favored by the consumers. The crossbreeding technology has been widely applied in animal breeding, as well as these hybrid individuals have been farmed for a long time forming gradually new indigenous breeds. Therefore, it is highly significant to understand the mechanisms regulating muscle growth and development of the generated crossbreed. The longissimus dorsi (LD) muscle is one of the important representatives of muscle tissues in animals, which is closely related to the individual skeletal muscle growth and development, as well as with the intramuscular fat content, and tenderness. A series of reported myogenic regulatory molecules regulate the myoblast’s proliferation and myophagism, further influencing the growth process, including the identified genes MYOG (Rudnicki and Jaenisch, 1995) and MYF5 (Xu et al., 2019), which are involved in regulating the myoblast differentiation, muscle growth, and meat quality traits in livestock. Previous studies have indicated that MYOD as a key regulator of myotube formation promotes the myotube’s differentiation (Tapscott, 2005; Wang et al., 2017). The CAPN family are important candidate genes in the growth and degradation of the muscle fibers, and are specifically expressed in the skeletal muscle (Gandolfi et al., 2011).
In fact, the studies indicated that the biological processes were not only regulated by protein-coding RNA—mRNA, as the sequencing technology developed. Non-coding RNAs (ncRNAs), including long non-coding RNA (lncRNA), circular RNA (circRNA), and microRNA (miRNA), are profoundly involved in diverse biological processes and are regulating them by various mechanisms. The lncRNAs universally acknowledged participating in chromatin transcriptional/epigenetic regulation by interacting with the chromatin regulators as “molecular scaffold” or decoys to activate or repress transcription (Caretti et al., 2006; Korostowski et al., 2012). Many lncRNAs have been proven to play a vital role in skeletal muscle development; for example, the lncRNA MAR1 positively correlates with muscle differentiation and growth in vitro and in vivo (Zhang et al., 2018). The lnc-smad7/miR-125b/Smad7 (SMAD family member 7) and IGF2 axes are instrumental in myoblast differentiation and regeneration of muscle in two different pig breeds (Song et al., 2018). In addition, continuously growing discoveries have reported the role of novel circRNA in skeletal muscle. CircZfp609 derived from Zinc Finger Protein 609, can inhibit the myogenic differentiation via the sponge miR-194-5p in the mouse myoblast cell line C2C12 (Wang et al., 2019b). The circFGFR4 was generated from the fibroblast growth factor receptor 4 (FGFR4) and could simulate the bovine primary myoblast’s differentiation through the circFGFR4-miR-107-WNT3A axis in cattle (Li et al., 2018). MicroRNA response elements are considered to be “talking mediators” of mRNAs, lncRNAs, and transcribed pseudogenes (Ji et al., 2019), and these response elements have important roles in various biological processes by forming a large number of complex regulatory networks. Numerous reports on the conjoint effect of miRNA and mRNA conjointly on developing skeletal muscle have been proved. FGFR1, which could prevent muscle fibrogenesis, is a functional target of miR-214-3p (Arrighi et al., 2021). The miR-183 and miR-96 were found to negatively regulate fat usage in the skeletal muscle via targeting FoxO1 and PDK4 (Wang et al., 2021a). Zhang et al. (2021) experimentally confirmed that miR-22-3p regulated the WFIKKN2 gene in adipocyte differentiation in muscle fat metabolism of Yanbian cattle. A targeted relationship between the oar-miR-655-3p and oar-miR-381-5p with ACSM3 and ABAT has been found to have crucial roles in sheep muscle organogenesis, myoblast migration (Sun et al., 2019).
In recent years, with the deepening of the research on the function of miRNAs, a new theory named competing endogenous RNA (ceRNA) has emerged. At the same time, some studies reported that mRNAs, lncRNAs, and circRNAs might regulate the gene function via miRNA and act as ceRNAs in various biological processes (Salmena et al., 2011; Yu et al., 2019). In the whole transcriptome, a comprehensive post-transcriptional regulatory network formed by ceRNA activity has greatly widened the cognition of functional genetic information in the genome. There are effective interactions among the lncRNA, circRNA, and mRNA with miRNA, and they can take significant effect in various processes of regulation in animals. LncRNAs act as molecular sponges for the miRNAs that specifically inhibit the target mRNAs so that they can give play to the protection of mRNAs (Li et al., 2019). For instance, a previous study reported that MAML1 and MEF2C, as transcription factors, activate the late-differentiation muscle genes, and linc-MD1 can regulate their expression as a ceRNA by sponging miR-133 and miR-135 (Cesana et al., 2011). LncRNA H19 can regulate muscle differentiation as a molecular sponge for the LET7 family in the developing embryo and adult muscles (Kallen et al., 2013).
Until now, most studies have been based on focusing on the cattle-yak for exploring the male sterility mechanism and barely referred to the superiority in the growth mechanism. Here, we have measured the growth traits of cattle-yaks and domestic Ashidan yaks under the same feeding and management, and systematically explored the differences of the longissimus dorsi muscles for the first time using the whole-transcriptome sequencing. Furthermore, the ceRNA network was constructed to identify the key factors involved in muscle growth and development. This study will thus help in improving yak breeding and provide new ideas for studying the genetic mechanism of muscle growth.
Materials and MethodsEthics Approval
All the animal experiments were approved by Lanzhou Institute of Husbandry and Pharmaceutical Sciences of the Chinese Academy of Agricultural Sciences (CAAS) with the grant number: No. 2019-002. All the slaughter as well as sampling procedures strictly complied with the Guidelines on the Ethical Treatment of Experimental Animals of China.
Phenotypic Data Collection and Samples Preparation
Thirty cattle-yaks (Aberdeen Angus ♂ × Yak ♀) and 30 yaks were tracked to measure the production performance indices (withers height, body weight, chest girth, and body length) at three stages of growth (i.e., birth, 3 months, and 6 months). Cattle-yaks (n = 3, 6 months old) and yaks (n = 3, 6 months old) were selected randomly to be slaughtered for longissimus dorsi muscle. These samples were collected for transcriptome sequencing and Real-time quantitative polymerase chain reaction analysis (RT-qPCR). All the samples were stored in liquid nitrogen (−80°C) for the subsequent tests.
RNA Isolation and Illumina Sequencing
The total RNA was isolated with TRIzol (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s instructions, and the concentration and quality of RNA were assessed by 1.5% agarose gel electrophoresis and Thermo Scientific NanoDrop 2000c (ThermoFisher Scientific Inc., Waltham, MA, United States).
Equal quantities of RNA were pooled from each sample. Then, the TruSeq Stranded Total RNA with Ribo-Zero Gold Kit (Illumina, San Diego, CA, United States) was used for digesting the ribosomal RNA (rRNA) in the DNA-free RNA. According to the manufacturer’s instructions, we performed the construction of library preparation with NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ipswich, MA, United States). The size and purity of libraries were validated by Agilent Technologies 2100 Bioanalyzer (Agilent, Santa Clara, CA). Finally, the samples were sequenced using Illumina HiSeq 2500 Technology (LC Sciences, Houston, TX, United States) with a 150-bp paired-end run.
Data Preprocessing, Read Mapping, and Transcript Assembly
Raw reads generated during high-throughput sequencing were in fastq format. The raw sequencing dataset supporting the results of cattle-yaks in this study was deposited at NCBI’s Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/, accession number: PRJNA753699). The raw sequencing dataset of yaks has been uploaded to the NCBI Database in the previous study (Ma et al., 2020). To get high-quality clean reads that could be used for later analysis, the Trimmomatic software (Bolger et al., 2014) was applied to remove adaptors, low-quality bases, and N-bases.
The remaining high-quality cleaned reads of each sample were aligned to the yak reference genome (BosGru_v2.0) using the HiSAT2 software (Kim et al., 2015). All the samples were assessed by genomic and gene alignment. After obtaining the comparison result bam file, the reads on the comparison gene were assembled using the StringTie software (Pertea et al., 2015), and every single transcript assembled by each sample was fused and spliced into a merged transcript.
LncRNAs Identification and LncRNAs Target Gene Prediction
We performed the following steps to obtain the potential lncRNA candidates for subsequent analysis: (1) The Cuffcompare software (Ghosh and Chan, 2016) was used to compare the merged transcripts with the reference transcripts one by one, and the transcripts marked with “i,” “u,” “x,” and “o” were retained after clarifying the position class of the remaining transcripts. (2) Transcripts with length >200 nt and exon number ≥2 were obtained. (3) The above transcripts were analyzed about the coding ability by CPC2 (Kang et al., 2017), CNCI (Sun et al., 2013), PLEK (Li et al., 2014), and Pfam (Sonnhammer et al., 1998) software to remove the transcripts with coding potential and obtain potential lncRNA candidates.
Cis-acting and trans-acting modes are two main ways to predict the targets of lncRNAs (Mercer et al., 2009). We calculated the locations of the paired lncRNAs and mRNAs for the cis-acting prediction. The lncRNA with no nearest protein-coding gene within 100 kb upstream or downstream was excluded in subsequent analysis. For trans-acting regulatory mode, the LncTar was used to calculate the free energy between them to predict the regulatory targets, since the expression of lncRNA is determined to be independent of the location of mRNA.
CircRNA Identification
We used the CIRI software and predicted the circRNA based on the BWA software (Li and Durbin, 2009; Gao et al., 2015). Since it is an authoritative software, it has the characteristics of high sensitivity and multiple screening for reducing false positives. Firstly, we aligned the clean reads to the reference genome to obtain the SAM file using the BWA software. Then, the CIRI software was used to scan for PCC signals (paired chiastic clipping signals), and circRNA sequences were predicted based on junction reads and GT-AG cleavage signals. The expressional levels of circRNAs were quantified by the RPM algorithm.
Differential Expression Genes and Pathway Analysis
The expression levels of the mRNAs and lncRNAs were calculated through the fragments per kilobase of transcript per million reads mapped (FPKM) value (Liu et al., 2018) using the Cuffdiff program. The DESeq software (Anders and Huber, 2010) was used to standardize the counts of each sample and calculate the fold change (FC). A negative binomial distribution test (NB) was used to test the difference significance of counts. The differentially expressed mRNAs (DEMs), lncRNAs (DELs), and circRNAs (DECs) were screened according to the results of |log 2 FC| >1 and p < 0.05 eventually.
To further understand gene function, Gene Ontology (GO, http://www.geneontology.org) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/) were used for enrichment analysis of functional pathways. Each GO and KEGG enrichment term was confirmed by Hypergeometric Distribution Test. Then, the p-value was corrected by Benjamini and Hochberg multiple tests. The enrichment with p-value lower than 0.05 was considered significant.
Construction of ceRNA Network
To acquire a better understanding of the interactions of the mRNAs, lncRNAs, circRNAs, and miRNAs, a lncRNA–circRNA–miRNA–mRNA regulatory network was constructed based on the ceRNA hypothesis (Salmena et al., 2011). MiRanda (John et al., 2004) was used to predict the pairs of miRNA–lncRNA, miRNA–mRNA, and miRNA–circRNA. The Spearman correlation coefficient (SCC) was used to evaluate the pairwise correlations of miRNA–lncRNA, miRNA–mRNA, and miRNA–circRNA; the value greater than 0.8 was considered relevant for constructing the network and p < 0.05 was regarded as being statistically significant. The Cytoscape software (version 3.5.1) was used to display the results visually.
Quantitative Real-Time PCR for Validating Gene Expression
The RNA for verifying the gene expression was the same as RNA used in the above Illumina sequencing. The RNA was reverse transcribed into cDNA by HiScript® II 1st Strand cDNA Synthesis Kit (Vazyme, Nanjing, Jiangsu, China). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), as endogenous control, was used to normalize target gene expression. The total 20-μl reaction system of qPCR included 10 ng of cDNA, 10 μl of SYBR Premix Ex Taq II (TaKaRa, Dalian, China), each 1 μl of forward primer and reverse primer, and 7 μl of ddH2O. The RT-qPCR conditions included the preincubation at 95°C for the 30 s, 45 cycles of 10 s at 95°C and 60 s at 59°C, then ended at 72°C for 30 s. Each experiment was repeated three times, and the results of relative RNA expression were calculated according to the cycle threshold (Ct) value 2–ΔΔCt (Livak and Schmittgen, 2001).
ResultsComparison of Production Performance
Analysis of the measurement results showed that the cattle-yaks had a significant improvement in production performance compared to the yaks in each age stage (p < 0.001). As can be seen from Table 1, cattle-yaks were stronger and taller at birth than yaks. In terms of body weight, at birth, the cattle-yak increased by about 22.60% compared to the yaks, and the body length, height, and chest girth of cattle-yaks were also increased by 17.62%, 11.06%, and 8.92%, respectively. From birth to 6 months old, the cattle-yak showed obvious heterosis with varying degrees of improvement in various indicators.
The production performance between cattle-yaks and yaks.
Items
Birth (Mean ± SE)
3 months (Mean ± SE)
6 months (Mean ± SE)
Yaks
Cattle-yaks
p-Value
Yaks
Cattle-yaks
p-Value
Yaks
Cattle-yaks
p-Value
Body weight (kg)
16.37 ± 0.23
20.70 ± 0.41
<0.001
37.56 ± 0.51
48.21 ± 0.62
<0.001
81.40 ± 1.31
90.76 ± 0.81
<0.001
Body height (cm)
57.23 ± 0.87
63.56 ± 0.56
<0.001
73.73 ± 0.98
84.80 ± 0.95
<0.001
85.43 ± 1.12
91.73 ± 0.86
<0.001
Body length (cm)
50.50 ± 1.10
59.40 ± 0.55
<0.001
75.63 ± 0.72
87.27 ± 0.78
<0.001
89.90 ± 0.77
98.07 ± 1.04
<0.001
Chest girth (cm)
58.63 ± 0.79
63.86 ± 0.54
<0.001
84.46 ± 0.60
97.23 ± 0.68
<0.001
112.83 ± 1.23
120.77 ± 1.20
<0.001
Note: SE: standard error.
The Prediction and Characteristic of LncRNAs and CircRNAs
After filtering out the low-quality and supernumerary reads, we obtained a total of 310,089,232 and 308,742,564 clean reads with greater than 93.77% of Q30 from the LD muscles of cattle-yaks and yaks, respectively. As shown in Supplementary Table S1, the successful alignment of reads to the reference genome was approximately over 94.55%.
After rigorous screening and filtration, we detected a large number of lncRNAs and cirRNAs (Figures 1A,E). Among them, 1,817 lncRNAs with an average length of 1,449 bp were discovered as novel lncRNAs (Figure 1B). The largest proportion of lncRNAs was over 2,000 bp (19.37%) and the proportion of lncRNA containing two exons was about 70.61% (Figure 1C). The total lncRNAs of 135 exonic antisense, 219 intronic antisense, 105 intergenic downstream antisense, 202 intergenic upstream, 204 exonic sense, 380 intronic sense, 141 intergenic downstream sense and 140 intergenic up-stream sense were identified in our results (Figure 1D). The average length of detected circRNAs was 3,250 bp and the length of the most circRNAs was over 2,000 bp (16.92%) (Figure 1F). As evident in Figure 1E, the circRNAs of sense-overlapping accounts for 89% of the total. The number of circRNAs located in the exonic and intronic was 332 and 137, respectively, while the antisense-overlapping circRNAs occupied about 1%.
The description of identified lncRNAs and circRNAs. (A) Screening of the candidate lncRNAs in longissimus dorsi muscle. (B) The length distribution of the novel lncRNAs. (C) Exon number distribution of novel lncRNAs. (D) The classification of novel lncRNAs. (E) The structure type pie chart of circRNAs. (F) The length distribution of circRNAs.
Differentially Expressed mRNAs, lncRNAs, and circRNAs Between CY and Y Groups
The study identified 7,126 mRNAs, 791 lncRNAs, and 1,057 circRNAs to have significant differential expressions (Supplementary Tables S2–S4). There were 6902 DE mRNAs, 742 DE lncRNAs, and 273 DE circRNAs, respectively, which were detected in the cattle-yaks and yaks (Figures 2A–C). However, 119 mRNAs, 41 lncRNAs, and 232 circRNAs were detected to express only in the CY group. Other 105 mRNAs, 8 lncRNAs, and 552 circRNAs were identified only in the Y group. To find the overall distribution of differential expression, the volcano plot was drawn based on the results of the differential expression. Compared with the yaks of the control group, we found 3,563 upregulated mRNAs, 455 upregulated lncRNAs, and 353 upregulated circRNAs in the cattle-yaks. Meanwhile, there were 3,563 downregulated mRNAs, 336 downregulated lncRNAs, and 704 downregulated circRNAs, respectively, in the cattle-yaks (Figures 2D–F).
Comparative analysis of the differentially expressed mRNAs and ncRNAs between cattle-yaks and yaks. The specific and shared mRNAs (A), lncRNAs (B), and circRNAs (C) between the two groups. The analysis of differentially expressed mRNAs (D), lncRNAs (E), and circRNAs (F) between two groups. Note: The red and green points represented upregulated and downregulated mRNAs, lncRNAs, and circRNAs, respectively. The gray points represented no significant differences. The gray vertical lines showed |log2FC| = 1, and the gray horizontal lines showed p = 0.05.
Functional Analysis
GO analysis described the molecular functions performed by the DE ncRNAs, the cellular environment in which they were located, and the biological processes involved. We respectively selected the top 10 biological processes with the most significant enrichment of mRNAs, lncRNAs, and circRNAs (Figure 3A). Obviously, the DEMs were mainly enriched in regulation of osteoblast proliferation, negative regulation of fat cell proliferation, and regulation of cellular response to hypoxia. The targeted genes of DELs were most enriched in terms involving regulation of fat cell differentiation, ephrin receptor signaling pathway, and muscle contraction. Protein autophosphorylation, fatty acid catabolic process, and regulation of Rho protein signal transduction were the most significant enrichment terms for the sourced genes for DECs. As evident from Figures 3B–D, we conducted the KEGG analysis using KEGG public pathway database and draw the augmented scatter diagram of the selected target genes. Between cattle-yaks and yaks, the DEMs were enriched in the PI3K−Akt signaling pathway, MAPK signaling pathway, Fatty acid metabolism, Citrate cycle, etc. (Figure 3B). The DEL adjacent genes were significantly related to the protein digestion and absorption, GnRH signaling pathway, alanine, aspartate and glutamate metabolism, and so on (Figure 3C). The host genes of DECs were significantly associated with some pathways, such as those related to vitamin digestion, ABC transporters, cGMP−PKG signaling pathway, etc. (Figure 3D). Interestingly, some DEMs and the host genes of DECs were found to be enriched in the same pathways (i.e., hippo signaling pathway, cell cycle-caulobacter, and calcium signaling pathway), and some DEM and DEL adjacent genes were enriched in the same oxytocin signaling pathways.
(A) GO analysis with the top 10 enrichment biological processes for DEMs, DELs targets, and DECs host genes between cattle-yaks and yaks. KEGG analysis with the top 20 KEGG enrichment pathways for DEMs (B), DELs target genes (C), and DECs host genes (D) between the cattle-yaks and yaks.
In order to further explore the DEMs related to muscle growth and fatness, we screened out the related genes. These 117 DEMs related to muscle development and fat deposition are shown in Supplementary Table S5. The functional predictions of GO in DEMs (Figure 4A) mainly focused on some terms of skeletal muscle development, muscle cell differentiation, and regulation of canonical Wnt signaling pathway, including positive regulation of myoblast differentiation and skeletal muscle fiber development. According to the KEGG pathway analysis (Figure 4B), some pathways were enriched significantly by these DEMs, such as Hippo signaling pathway, regulating many biological processes involved in proliferation, survival and differentiation of cell, organ size, and tissue homeostasis. In addition, some DEMs enriched significantly in Notch signaling pathway that can regulate the differentiation and development of cells, tissues, and organs through the interaction between the adjacent cells.
Functional analysis of 117 differentially expressed genes associated with muscle development and fatness between the cattle-yak and yak. GO (A) and KEGG (B) analysis.
Construction of the ceRNA Coregulatory Network
Previous studies have shown that mRNAs, lncRNAs, and circRNAs may regulate gene function through miRNAs as ceRNAs in different processes (Salmena et al., 2011; Yu et al., 2019), indicating that ceRNAs and their miRNAs may work in concert with each other. Combined with the DEMs, DELs, and DECs related to muscle development and co-differentially expressed, we constructed the integrated ceRNA network. This ceRNA network contained 11 DEMs, 6 DELs, 8 DECs, and 33 relationships (Figure 5). Tcons-00034903 and bta-miR-2039 have a shared target gene WNT4. Similarly, we also found the same results in bta-miR-2316-Tcons-00029868-MYH4 and bta-miR-1777a-Tcons-00027748-SIX5. This ceRNA network might provide valuable information for the development of the longissimus dorsi in cattle-yaks and yaks.
The ceRNA co-regulation network. The blue circle, yellow box, red triangle and green V-type represented the differentially expressed mRNAs, lncRNAs, circRNAs, and miRNA, respectively. The solid line indicated the co-regulation between miRNAs and other transcripts. The dotted line indicated the co-regulation between the lncRNAs and mRNAs.
Real-Time Quantitative PCR Validation of Sequencing Data
Four mRNAs (MYH14, LOC106700760, PIK3R2, and FGFR4), 2 lncRNAs (Tcons-00034903 and Tcons-00004303), and 2 circRNAs (circ00012096 and circ00012564) were selected randomly for verifying the sequencing results through real-time quantitative PCR (Figure 6). The primers were designed using the Primer-BLAST web tool from the National Center for Biotechnology Information (NCBI) (Supplementary Table S6). Their expression patterns were highly consistent with the sequencing results, indicating that the gene expression profiles obtained in this study had high repeatability and reliability.
The results of the real-time quantitative PCR validation of the expression level. ** indicates p < 0.01, * indicates p < 0.05. The data represented the mean ± SEM from three biological replicates, and each measurement was repeated 3 times.
Discussion
Muscle growth is a complex economic trait owing to various physiological and biochemical indices. It concomitantly involves many gene expressions and regulations. To date, people have a great demand for high-quality meat and nutrition with improving living standards; it is an urgent problem for improving the performance of yaks. Cattle-yak, as hybrid offspring of yak (♀) and cattle (♂), has been found to have significant improvement in the production performance (Guo et al., 2019; Dingkao et al., 2020; Luo et al., 2020; Jiang et al., 2021). However, the economic benefit of cattle-yaks has been ignored for a long time. The variations in the genes and proteins in the muscle structure are presumed generally to be affected by production performance (Joo et al., 1999); this study is the first time to systematically explore the differences of growth and development between cattle-yaks and yaks based on transcriptomics. To study the regulatory mechanism of growth and muscle development of cattle-yaks and yaks in a better way, we performed transcriptome analysis of the longissimus dorsi muscle. It was also the first time to compare the differences in the expression profiles of mRNA, lncRNA, and circRNA between the cattle-yaks and yaks, so as to determine the key factors involved in muscle growth and development. We measured the phenotypic data strictly of 30 cattle-yaks and yaks at three age groups of growth using the standard method of measurement, and the results obtained were similar to the previous studies; in each period of cattle-yaks with the same feeding conditions, production performance indexes of cattle-yaks were significantly higher than those of yaks (p < 0.001). We speculated that these differences may be due to genetic factors rather than the effects of feeding management on these traits. Therefore, the study on the mechanism regulating muscle growth and development in hybrid yaks and yak breeds can better help in enhancing the production performance of the yaks, providing help for exploring the regulatory mechanism of muscle development in mammals.
Typically, the growth and development of muscles are regulated by the core genes and signal transduction pathways (Myers et al., 2006; Ayuso et al., 2015). Compared to the yak group, a total of 7,126 DEMs, 791 DELs, and 1,057 DECs were identified in the cattle-yak group. Subsequently, the GO and KEGG pathway enrichment analysis revealed some important DEMs related to muscle growth and fat deposition, which was consistent with the results obtained after measuring the production performance, indicating that the production performance of the cattle-yak was better than that of the yak (Tumennasan et al., 1997). In addition, some DEMs, DELs, and DECs related to the immune system have reflected the adaptation of the yak to the high-altitude environment, and the adaptation of the yaks to the cold and high-altitude environment was well-preserved by the cattle-yaks. In succession, we focused on some DEMs related to the production performance of cattle-yaks and yaks. A total of 117 DEMs were identified that were related to the differentiation and proliferation of myoblasts, AMPK signaling pathway, MAPK signaling pathway, PI3K–Akt signaling pathway, etc. (Neri et al., 2002; Anderson, 2006; Gehart et al., 2010; Thomson, 2018). Among these DEMs, some have been found to have known functions in muscle growth and development. For example, myostatin (MSTN), as a member of the TGF-β superfamily, has a proven role as a growth differentiation factor (Kollias and McDermott, 2008; Li et al., 2008) playing an important role in mice (McPherron et al., 1997), cattle (McPherron and Lee, 1997), and humans (Schuelke et al., 2004) by negatively regulating the growth and differentiation of myoblasts, as well as other mammals via controlling both the activation and proliferation of the satellite cells (skeletal muscle stem cells) (McCroskery et al., 2003). The expression level of the MSTN gene in our study was found to be downregulated in the cattle-yaks; hence, the expression of this gene was considered to be consistent with its function. This also provided new evidence for the function of the MSTN gene in longissimus dorsi of cattle-yaks and the high conservation across species. Myogenin (MYOG), including myogenic factor 6 (MYF6), is a regulatory factor in the family of Myogenic regulatory factors (MRFs), which is mainly involved in the fusion and differentiation of myoblasts (CHARGÉ and RUDNICKI, 2004; Buckingham and Vincent, 2009). Interestingly, the paired box transcription factor 7 (PAX7) gene and the MYOG gene were found to be downregulated, but the MYF6 expression was upregulated in the cattle-yaks. Previous studies have shown that overexpression of gene PAX7 can induce MYOG expression to inhibit myogenesis and prevent the differentiation of the muscle cell (Chen et al., 2010; Dey et al., 2011). Therefore, this also justified the downregulation of both PAX7 and MYOG gene expression in cattle-yaks with better production performance than yaks. Additionally, some genes, such as actin alpha 1 skeletal muscle (ACTA1) and actin alpha cardiac muscle 1 (ACTC1), also play a key role in the differentiation and fusion of muscle cells, and have a positive impact on the myogenesis of the skeletal muscles (Ciecierska et al., 2020). Although these studies have helped in predicting the functions and accuracy of the key genes, the other functional DEMs related to muscle growth and development still need to be further studied on their expression regulation in longissimus dorsi muscle.
The main non-coding RNAs, the lncRNAs and circRNAs, are receiving more and more attention, as they can participate in the regulation of various biological processes in different ways (Liu et al., 2017; Marchese et al., 2017; Wang et al., 2019a). The identified DELs and DECs in this study were involved in regulating the promotion of muscle growth and development. LncRNA can exert its important action through various biological and pathological processes by demonstrating trans, cis, and antisense effects. Our results showed abundant differentially expressed lncRNAs in the skeletal muscle of cattle-yak and yak, suggesting that the lncRNAs may not be the exclusive by-products of mRNA in the cattle-yak but have specific roles. The long non-coding RNA acts differently in the nucleus and cytoplasm due to their different locations, where they are involved in the muscle development regulation on both embryonic and growing stages. According to the KEGG analysis, the differentially expressed lncRNAs and circRNAs were identified to functionally relate to some hormone regulation, myoblast proliferation, and metabolic pathways. CircRNAs being another type of non-coding RNA also act as an important regulatory role in skeletal muscle growth and development (Greco et al., 2018). Multiple studies have confirmed that abundant circRNAs exist in the skeletal muscle and their expression levels change dynamically during the process of myoblast differentiation (Legnini et al., 2017). Notably, we found that many DECs host genes are enriched in the HIF-1 signaling pathway in the KEGG analysis. HIF-1 signaling pathway is the core signaling pathway induced by hypoxia involved in regulating the proliferation and differentiation of myoblasts upon hypoxic conditions (Ogilvie et al., 2000). HIF-1 pathway can upregulate its target genes, some of which involve the regulation of proliferation and differentiation of myoblasts. These DEC host genes were found to be significantly enriched in the HIF-1 signaling pathway, indicating a close relationship with the regulation of muscle growth and development under high-altitude and hypoxic environments. The regulatory mechanism of these circRNAs and their host genes on muscle development hence deserves further studies. Furthermore, our studies also indicated some lncRNAs and circRNAs to be specifically or mainly expressed in the longissimus dorsi muscle of cattle-yaks (e.g., Tcons-00005361, Tcons-00037918, and circRNA-02529), indicating that these ncRNAs were generated on purpose to have specific effects in the muscle development of the cattle-yaks.
In recent years, extensive studies on the function of miRNAs have provided a new theory named competing for endogenous RNA (ceRNA). To understand the process of development of the skeletal muscle, a ceRNA regulatory network was constructed based on the combination with the DEMs, DELs, and DECs related to muscle development and were expressed co-differentially. The results showed that 11 DEMs, 6 DELs, and 8 DECs cross-talked with another through 8 differential expressions of the microRNAs. This also indicated that the development of the longissimus dorsi muscle in cattle-yak was a complex regulation process of a balanced level of gene expression under a high-altitude and hypoxic environment. As reported, peroxisome proliferator-activated receptor delta (PPARD) acted as a vital regulator in adipogenesis and lipid metabolism (Kim et al., 2006; Chen and Yang, 2014). Ankyrin repeat domain 6 (ANKRD6) also played a role in regulating crucial events in developing vertebrates and invertebrates (Schwarz-Romond et al., 2002; Moeller et al., 2006). Therefore, we speculated that these ncRNAs might also contribute to muscle development by indirectly regulating the gene expression of PPARD and ANKRD6. Not only that, we observed three important ceRNA subnetworks from the ceRNA network, showing that TCONS-00024051 and its target SOX8 “talked” to each other through the same bta-miR-1777a response element, while TCONS-00034903 and its target WNT4, TCONS-00029868 and its target ACTA1 “talked” to each other through bta-miR-12039 and bta-miR-4449 response elements, respectively. Therefore, we speculated that these three subnetworks may be crucially associated and function in regulating muscle development. Notably, as a member of the “unconventional” non-muscle myosin II family of molecular motors, myosin heavy chain 14 (MYH14) gene has been identified as a key regulator of muscle fiber type (van Rooij et al., 2009; Bell et al., 2010). In our results, MYH14 and Tcons-00029868 were both upregulated in longissimus dorsi muscle, which indicates that this lncRNA may have a cis-regulatory relationship with MYH14 gene. From our results, these DEMs, DELs, and DECs were not only involved in the process of muscle development, but might also be involved in lipogenesis as ceRNAs.
Conclusion
In conclusion, we compared the production performance of cattle-yaks and yaks, it was also the first time to compare the expressional features of mRNAs, lncRNAs, and circRNAs in the longissimus dorsi tissue of cattle-yaks and yaks. In our results, the abundant mRNAs and ncRNAs were identified, some of which were specifically expressed in cattle-yaks. According to the bioinformatics analyses, the ncRNAs were found to be not only connected with the myoblast differentiation and proliferation, skeletal development, and signaling pathway of muscle growth, but be also useful as ceRNAs of important transcription factors (such as SOX8 and PPARD). In addition, the ceRNA network (11 DEMs, 6 DELs, and 8 DECs) that we constructed may have the considerable effects on regulation of muscle growth. This study provided new insights into the genetic basis of muscle growth and laid the foundation for further study of the role of these ncRNAs in regulating muscle growth.
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 below: https://www.ncbi.nlm.nih.gov/, PRJNA753699.
Ethics Statement
The animal study was reviewed and approved by the Lanzhou Institute of Husbandry and Pharmaceutical Sciences Chinese Academy of Agricultural Sciences. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
CL, PY, and CH contributed to conception and design of the study. CH and FG organized the database. FG and XM performed the statistical analysis. CH, FG, and XM wrote the first draft of the manuscript. RFD, RQD, ZZ, and GB collected the production data. XG, PB, XM, XW, and MC wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.
Funding
Work presented in this paper was supported by the Project of Innovative Research on Yak Molecular Breeding Technology under the Innovation Population of Basic Research in Gansu Province, grant number 20JR5RA580, the Agricultural Science and Technology Innovation Program, grant number CAAS-ASTIP-2014-LIHPS-01, the National Beef Cattle Industry Technology and System, grant number CARS-37, and the Herbivorous Livestock Industry and Technical System of Gansu Province, grant number GARS-08.
Conflict of Interest
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.2021.772557/full#supplementary-material
ReferencesAndersS.HuberW. (2010). Differential Expression Analysis for Sequence Count Data. Genome Biol.11 (10), R106. 10.1186/gb-2010-11-10-r106AndersonD. H. (2006). Role of Lipids in the MAPK Signaling Pathway. Prog. Lipid Res.45 (2), 102–119. 10.1016/j.plipres.2005.12.003ArrighiN.MoratalC.SavaryG.FassyJ.NottetN.PonsN. (2021). The FibromiR miR-214-3p Is Upregulated in Duchenne Muscular Dystrophy and Promotes Differentiation of Human Fibro-Adipogenic Muscle Progenitors. Cells10 (7), 1832. 10.3390/cells10071832AyusoM.FernándezA.NúñezY.BenítezR.IsabelB.BarragánC. (2015). Comparative Analysis of Muscle Transcriptome between Pig Genotypes Identifies Genes and Regulatory Mechanisms Associated to Growth, Fatness and Metabolism. PLoS One10 (12), e0145162. 10.1371/journal.pone.0145162BellM. L.BuvoliM.LeinwandL. A. (2010). Uncoupling of Expression of an Intronic microRNA and its Myosin Host Gene by Exon Skipping. Mol. Cel Biol.30 (8), 1937–1945. 10.1128/mcb.01370-09BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a Flexible Trimmer for Illumina Sequence Data. Bioinformatics30 (15), 2114–2120. 10.1093/bioinformatics/btu170BuckinghamM.VincentS. D. (2009). Distinct and Dynamic Myogenic Populations in the Vertebrate Embryo. Curr. Opin. Genet. Dev.19 (5), 444–453. 10.1016/j.gde.2009.08.001CarettiG.SchiltzR. L.DilworthF. J.Di PadovaM.ZhaoP.OgryzkoV. (2006). The RNA Helicases P68/p72 and the Noncoding RNA SRA Are Coregulators of MyoD and Skeletal Muscle Differentiation. Dev. Cel11 (4), 547–560. 10.1016/j.devcel.2006.08.003CesanaM.CacchiarelliD.LegniniI.SantiniT.SthandierO.ChinappiM. (2011). A Long Noncoding RNA Controls Muscle Differentiation by Functioning as a Competing Endogenous RNA. Cell147 (2), 358–369. 10.1016/j.cell.2011.09.028ChargéS. B. P.RudnickiM. A. (2004). Cellular and Molecular Regulation of Muscle Regeneration. Physiol. Rev.84 (1), 209–238. 10.1152/physrev.00019.2003ChenJ.-F.TaoY.LiJ.DengZ.YanZ.XiaoX. (2010). microRNA-1 and microRNA-206 Regulate Skeletal Muscle Satellite Cell Proliferation and Differentiation by Repressing Pax7. J. Cel Biol.190 (5), 867–879. 10.1083/jcb.200911036ChenL.YangG. (2014). PPARs Integrate the Mammalian Clock and Energy Metabolism. PPAR Res.2014, 1–6. 10.1155/2014/653017CiecierskaA.MotylT.SadkowskiT. (2020). Transcriptomic Profile of Primary Culture of Skeletal Muscle Cells Isolated from Semitendinosus Muscle of Beef and Dairy Bulls. Ijms21 (13), 4794. 10.3390/ijms21134794DeyB. K.GaganJ.DuttaA. (2011). miR-206 and -486 Induce Myoblast Differentiation by Downregulating Pax7. Mol. Cel Biol.31(1),203–214. 10.1128/MCB.01009-10DingkaoR. Q.ShiH. M.LiP. X.CairangN. R.MaD. L.NiuX. L. (2020). Application of Early Pregnancy Diagnostic Technique in Xianan Cattle Production. China Cattle Sci.46 (5), 11–14. GandolfiG.PomponioL.ErtbjergP.KarlssonA. H.Nanni CostaL.LametschR. (2011). Investigation on CAST, CAPN1 and CAPN3 Porcine Gene Polymorphisms and Expression in Relation to post-mortem Calpain Activity in Muscle and Meat Quality. Meat Sci.88 (4), 694–700. 10.1016/j.meatsci.2011.02.031GaoY.WangJ.ZhaoF. (2015). CIRI: an Efficient and Unbiased Algorithm for De Novo Circular RNA Identification. Genome Biol.16 (1), 4. 10.1186/s13059-014-0571-3GehartH.KumpfS.IttnerA.RicciR. (2010). MAPK Signalling in Cellular Metabolism: Stress or Wellness?EMBO Rep.11 (11), 834–840. 10.1038/embor.2010.160GhoshS.ChanC.-K. K. (2016). “Analysis of RNA-Seq Data Using TopHat and Cufflinks,” in Plant Bioinformatics: Methods and Protocols. Editor EdwardsD. (New York, NY: Springer New York), 339–361. 10.1007/978-1-4939-3167-5_18GrecoS.CardinaliB.FalconeG.MartelliF. (2018). Circular RNAs in Muscle Function and Disease. Ijms19 (11), 3454. 10.3390/ijms19113454GuoS.MaD.LiB.BaoZ.ZhangT.XuG. (2019). Determination of Growth and Development Indexes of Cattle-Yak in Gannan alpine Pasture Area. China Herbivore Sci.039 (002), 73–75. 10.3969/j.issn.2095-3887.2019.02.021JiQ.CaiG.LiuX.ZhangY.WangY.ZhouL. (2019). MALAT1 Regulates the Transcriptional and Translational Levels of Proto-Oncogene RUNX2 in Colorectal Cancer Metastasis. Cell Death Dis.10 (6), 378. 10.1038/s41419-019-1598-xJiangH.ZhangC. F.ZhangQ.SunG. M.CaoH. W.CiJiY. D. (2021). Study on Yak Economic Crossbreeding and Growth Performance of F1 Generation. Tibet J. Agric. Sci.43 (02), 22–24. JohnB.EnrightA. J.AravinA.TuschlT.SanderC.MarksD. S. (2004). Human MicroRNA Targets. Plos Biol.2 (11), e363. 10.1371/journal.pbio.0020363JooS. T.KauffmanR. G.KimB. C.ParkG. B. (1999). The Relationship of Sarcoplasmic and Myofibrillar Protein Solubility to Colour and Water-Holding Capacity in Porcine Longissimus Muscle. Meat Sci.52 (3), 291–297. 10.1016/S0309-1740(99)00005-4KallenA. N.ZhouX.-B.XuJ.QiaoC.MaJ.YanL. (2013). The Imprinted H19 lncRNA Antagonizes Let-7 microRNAs. Mol. Cel52 (1), 101–112. 10.1016/j.molcel.2013.08.027KangY.-J.YangD.-C.KongL.HouM.MengY.-Q.WeiL. (2017). CPC2: a Fast and Accurate Coding Potential Calculator Based on Sequence Intrinsic Features. Nucleic Acids Res.45 (W1), W12–W16. 10.1093/nar/gkx428KimD. J.BilityM. T.BillinA. N.WillsonT. M.GonzalezF. J.PetersJ. M. (2006). Pparβ/δ Selectively Induces Differentiation and Inhibits Cell Proliferation. Cell Death Differ13 (1), 53–60. 10.1038/sj.cdd.4401713KimD.LangmeadB.SalzbergS. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods12 (4), 357–360. 10.1038/nmeth.3317KolliasH. D.McDermottJ. C. (2008). Transforming Growth Factor-β and Myostatin Signaling in Skeletal Muscle. J. Appl. Physiol.104 (3), 579–587. 10.1152/japplphysiol.01091.2007KorostowskiL.SedlakN.EngelN. (2012). The Kcnq1ot1 Long Non-coding RNA Affects Chromatin Conformation and Expression of Kcnq1, but Does Not Regulate its Imprinting in the Developing Heart. Plos Genet.8 (9), e1002956. 10.1371/journal.pgen.1002956LegniniI.Di TimoteoG.RossiF.MorlandoM.BrigantiF.SthandierO. (2017). Circ-ZNF609 Is a Circular RNA that Can Be Translated and Functions in Myogenesis. Mol. Cel66 (1), 22–37. e29. 10.1016/j.molcel.2017.02.017LiA.ZhangJ.ZhouZ. (2014). PLEK: a Tool for Predicting Long Non-coding RNAs and Messenger RNAs Based on an Improved K-Mer Scheme. BMC Bioinformatics15 (1), 311. 10.1186/1471-2105-15-311LiH.DurbinR. (2009). Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics25 (14), 1754–1760. 10.1093/bioinformatics/btp324LiH.WeiX.YangJ.DongD.HaoD.HuangY. (2018). circFGFR4 Promotes Differentiation of Myoblasts via Binding miR-107 to Relieve its Inhibition of Wnt3a. Mol. Ther. - Nucleic Acids11, 272–283. 10.1016/j.omtn.2018.02.012LiZ. B.KolliasH. D.WagnerK. R. (2008). Myostatin Directly Regulates Skeletal Muscle Fibrosis. J. Biol. Chem.283 (28), 19371–19378. 10.1074/jbc.M802585200LiZ.CaiB.AbdallaB. A.ZhuX.ZhengM.HanP. (2019). LncIRS1 Controls Muscle Atrophy via Sponging miR‐15 Family to Activate IGF1‐PI3K/AKT Pathway. J. Cachexia, Sarcopenia Muscle10 (2), 391–410. 10.1002/jcsm.12374LiuJ.LiuT.WangX.HeA. (2017). Circles Reshaping the RNA World: from Waste to Treasure. Mol. Cancer16 (1), 58. 10.1186/s12943-017-0630-yLiuX.LiuK.ShanB.WeiS.LiD.HanH. (2018). A Genome-wide Landscape of mRNAs, lncRNAs, and circRNAs during Subcutaneous Adipogenesis in Pigs. J. Anim. Sci Biotechnol9 (1), 76. 10.1186/s40104-018-0292-7LivakK. J.SchmittgenT. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods25 (4), 402–408. 10.1006/meth.2001.1262LuoQ. B.BurenC. G. T.SiQ.ZhaxiZ. M.ZhangS.ChangL. (2020). Study on Crossbreeding Effect between Angus Beef Cattle and Yak. China Cattle Sci.46 (5), 8–10. MaX.FuD.ChuM.DingX.WuX.GuoX. (2020). Genome-Wide Analysis Reveals Changes in Polled Yak Long Non-coding RNAs in Skeletal Muscle Development. Front. Genet.11 (365). 10.3389/fgene.2020.00365MarcheseF. P.RaimondiI.HuarteM. (2017). The Multidimensional Mechanisms of Long Noncoding RNA Function. Genome Biol.18 (1), 206. 10.1186/s13059-017-1348-2McCroskeryS.ThomasM.MaxwellL.SharmaM.KambadurR. (2003). Myostatin Negatively Regulates Satellite Cell Activation and Self-Renewal. J. Cel Biol.162 (6), 1135–1147. 10.1083/jcb.200207056McPherronA. C.LawlerA. M.LeeS.-J. (1997). Regulation of Skeletal Muscle Mass in Mice by a New TGF-P Superfamily Member. Nature387 (6628), 83–90. 10.1038/387083a0McPherronA. C.LeeS.-J. (1997). Double Muscling in Cattle Due to Mutations in the Myostatin Gene. Proc. Natl. Acad. Sci.94 (23), 12457–12461. 10.1073/pnas.94.23.12457MercerT. R.DingerM. E.MattickJ. S. (2009). Long Non-coding RNAs: Insights into Functions. Nat. Rev. Genet.10 (3), 155–159. 10.1038/nrg2521MoellerH.JennyA.SchaefferH.-J.Schwarz-RomondT.MlodzikM.HammerschmidtM. (2006). Diversin Regulates Heart Formation and Gastrulation Movements in Development. Proc. Natl. Acad. Sci.103 (43), 15900–15905. 10.1073/pnas.0603808103MyersS. A.WangS.-C. M.MuscatG. E. O. (2006). The Chicken Ovalbumin Upstream Promoter-Transcription Factors Modulate Genes and Pathways Involved in Skeletal Muscle Cell Metabolism. J. Biol. Chem.281 (34), 24149–24160. 10.1074/jbc.M601941200NeriL. M.BorgattiP.CapitaniS.MartelliA. M. (2002). The Nuclear Phosphoinositide 3-kinase/AKT Pathway: a New Second Messenger System. Biochim. Biophys. Acta1584 (2), 73–80. 10.1016/S1388-1981(02)00300-1OgilvieM.YuX.Nicolas-MetralV.PulidoS. M.LiuC.RueggU. T. (2000). Erythropoietin Stimulates Proliferation and Interferes with Differentiation of Myoblasts. J. Biol. Chem.275 (50), 39754–39761. 10.1074/jbc.M004999200PerteaM.PerteaG. M.AntonescuC. M.ChangT.-C.MendellJ. T.SalzbergS. L. (2015). StringTie Enables Improved Reconstruction of a Transcriptome from RNA-Seq Reads. Nat. Biotechnol.33 (3), 290–295. 10.1038/nbt.3122RudnickiM. A.JaenischR. (1995). The MyoD Family of Transcription Factors and Skeletal Myogenesis. Bioessays17 (3), 203–209. 10.1002/bies.950170306SalmenaL.PolisenoL.TayY.KatsL.PandolfiP. P. (2011). A ceRNA Hypothesis: The Rosetta Stone of a Hidden RNA Language?Cell146 (3), 353–358. 10.1016/j.cell.2011.07.014SchuelkeM.WagnerK. R.StolzL. E.HübnerC.RiebelT.KömenW. (2004). Myostatin Mutation Associated with Gross Muscle Hypertrophy in a Child. N. Engl. J. Med.350 (26), 2682–2688. 10.1056/NEJMoa040933Schwarz-RomondT.AsbrandC.BakkersJ.KühlM.SchaefferH.-J.HuelskenJ. (2002). The Ankyrin Repeat Protein Diversin Recruits Casein Kinase Iε to the β-catenin Degradation Complex and Acts in Both Canonical Wnt and Wnt/JNK Signaling. Genes Dev.16 (16), 2073–2084. 10.1101/gad.230402SongC.HuangY.YangZ.MaY.ChaogetuB.ZhuomaZ. (2019). RNA-seq Analysis Identifies Differentially Expressed Genes in Subcutaneous Adipose Tissue in Qaidaford Cattle, Cattle-Yak, and Angus Cattle. Animals9 (12), 1077. 10.3390/ani9121077SongC.WangJ.MaY.YangZ.DongD.LiH. (2018). Linc-smad7 Promotes Myoblast Differentiation and Muscle Regeneration via Sponging miR-125b. Epigenetics13 (6), 591–604. 10.1080/15592294.2018.1481705SonnhammerE.EddyS. R.BirneyE.BatemanA.DurbinR. (1998). Pfam: Multiple Sequence Alignments and HMM-Profiles of Protein Domains. Nucleic Acids Res.26 (1), 320–322. 10.1093/nar/26.1.320SunL.LuS.BaiM.XiangL.LiJ.JiaC. (2019). Integrative microRNA-mRNA Analysis of Muscle Tissues in Qianhua Mutton Merino and Small Tail Han Sheep Reveals Key Roles for Oar-miR-655-3p and Oar-miR-381-5p. DNA Cel Biol.38 (5), 423–435. 10.1089/dna.2018.4408SunL.LuoH.BuD.ZhaoG.YuK.ZhangC. (2013). Utilizing Sequence Intrinsic Composition to Classify Protein-Coding and Long Non-coding Transcripts. Nucleic Acids Res.41 (17), e166. 10.1093/nar/gkt646TapscottS. J. (2005). The Circuitry of a Master Switch: Myod and the Regulation of Skeletal Muscle Gene Transcription. Development132 (12), 2685–2695. 10.1242/dev.01874ThomsonD. (2018). The Role of AMPK in the Regulation of Skeletal Muscle Size, Hypertrophy, and Regeneration. Ijms19 (10), 3125. 10.3390/ijms19103125TumennasanK.TuyaT.HottaY.TakaseH.SpeedR. M.ChandleyA. C. (1997). Fertility Investigations in the F1 Hybrid and Backcross Progeny of Cattle (Bos taurus) Andyak (B. Grunniens) in Mongolia. Cytogenet. Cel Genet78 (1), 69–73. 10.1159/000134633van RooijE.QuiatD.JohnsonB. A.SutherlandL. B.QiX.RichardsonJ. A. (2009). A Family of microRNAs Encoded by Myosin Genes Governs Myosin Expression and Muscle Performance. Dev. Cel17 (5), 662–673. 10.1016/j.devcel.2009.10.013WangC.LiuW.NieY.QaherM.HortonH. E.YueF. (2017). Loss of MyoD Promotes Fate Transdifferentiation of Myoblasts into Brown Adipocytes. EBioMedicine16, 212–223. 10.1016/j.ebiom.2017.01.015WangH.MaM.LiY.LiuJ.SunC.LiuS. (2021a). miR‐183 and miR‐96 Orchestrate Both Glucose and Fat Utilization in Skeletal Muscle. EMBO Rep.22, e52247. 10.15252/embr.202052247WangJ.RenQ.HuaL.ChenJ.ZhangJ.BaiH. (2019a). Comprehensive Analysis of Differentially Expressed mRNA, lncRNA and circRNA and Their ceRNA Networks in the Longissimus Dorsi Muscle of Two Different Pig Breeds. Ijms20 (5), 1107. 10.3390/ijms20051107WangY.LiM.WangY.LiuJ.ZhangM.FangX. (2019b). A Zfp609 Circular RNA Regulates Myoblast Differentiation by Sponging miR-194-5p. Int. J. Biol. Macromolecules121, 1308–1313. 10.1016/j.ijbiomac.2018.09.039WangY.WangZ.HuR.PengQ.XueB.WangL. (2021b). Comparison of Carcass Characteristics and Meat Quality between Simmental Crossbred Cattle, Cattle‐yaks and Xuanhan Yellow Cattle. J. Sci. Food Agric.101 (9), 3927–3932. 10.1002/jsfa.11032XuX.MishraB.QinN.SunX.ZhangS.YangJ. (2019). Differential Transcriptome Analysis of Early Postnatal Developing Longissimus Dorsi Muscle from Two Pig Breeds Characterized in Divergent Myofiber Traits and Fatness. Anim. Biotechnol.30 (1), 63–74. 10.1080/10495398.2018.1437045YuC.LongfeiL.LongW.FengZ.ChenJ.ChaoL. (2019). LncRNA PVT1 Regulates VEGFC through Inhibiting miR‐128 in Bladder Cancer Cells. J. Cel Physiol.234 (2), 1346–1353. 10.1002/jcp.26929ZhangJ. S.XuH. Y.FangJ. C.YinB. Z.WangB. B.PangZ. (2021). Integrated microRNA-mRNA Analysis Reveals the Roles of microRNAs in the Muscle Fat Metabolism of Yanbian Cattle. Anim. Genet.52, 598–607. 10.1111/age.13126ZhangZ.-K.LiJ.GuanD.LiangC.ZhuoZ.LiuJ. (2018). A Newly Identified lncRNA MAR1 Acts as a miR-487b Sponge to Promote Skeletal Muscle Differentiation and Regeneration. J. Cachexia, Sarcopenia Muscle9 (3), 613–626. 10.1002/jcsm.12281