- Gansu Key Laboratory of Herbivorous Animal Biotechnology, Faculty of Animal Science and Technology, Gansu Agricultural University, Lanzhou, China
MicroRNAs (miRNAs) are a class of small non-coding RNAs that have been shown to play important post-transcriptional regulatory roles in the growth and development of skeletal muscle tissues. However, limited research into the effect of miRNAs on muscle development in goats has been reported. In this study, Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats with significant phenotype difference in meat production performance were selected and the difference in Longissimus dorsi muscle tissue expression profile of miRNAs between the two goat breeds was then compared using small RNA sequencing. A total of 1,623 miRNAs were identified in Longissimus dorsi muscle tissues of the two goat breeds, including 410 known caprine miRNAs, 928 known species-conserved miRNAs and 285 novel miRNAs. Of these, 1,142 were co-expressed in both breeds, while 230 and 251 miRNAs were only expressed in LC and ZB goats, respectively. Compared with ZB goats, 24 up-regulated miRNAs and 135 miRNAs down-regulated were screened in LC goats. A miRNA-mRNA interaction network showed that the differentially expressed miRNAs would target important functional genes associated with muscle development and intramuscular fat deposition. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed that the target genes of differentially expressed miRNAs were significantly enriched in Ras, Rap 1, FoxO, and Hippo signaling pathways. This study suggested that these differentially expressed miRNAs may be responsible for the phenotype differences in meat production performance between the two goat breeds, thereby providing an improved understanding of the roles of miRNAs in muscle tissue of goats.
Introduction
MicroRNAs (miRNAs) are a class of non-coding small RNA molecules (~22 nucleotides), which are evolutionarily conserved in eukaryotes (1). In recent years, miRNAs are well-recognized as negative regulators of gene expression at post-transcriptional level, in that they can either inhibit translation or promote degradation of mRNA by complementary binding to the 3′-untranslated regions (3′-UTR) of the target genes. The miRNAs are therefore involved in a wide variety of cell biological processes, including proliferation, differentiation, death, and fate specification (1, 2).
In modern animal husbandry, skeletal muscle is considered as the most economically important tissue of producing-meat livestock. Many studies have confirmed that miRNAs played essential roles in the growth and development of skeletal muscle. For example, miR-1 and miR-206 have been reported to facilitate differentiation and inhibit proliferation of bovine skeletal muscle satellite cells by targeting PAX7 (3). The over-expression of miR-486 induced skeletal muscle hypertrophy of mice via activation of protein kinase B (Akt) (4). The miR-27b regulated caprine myogenic proliferation and differentiation of skeletal muscle satellite cells by inhibiting the expression of PAX3 (5).
Up to now, research into expression profiles of miRNAs in the skeletal muscle tissue of domestic animals have mainly been focused on pigs (6–9), cattle (10–13), and sheep (14–17). It was found from these studies described above that miRNAs were differentially expressed in skeletal muscle at different developmental periods, or between different breeds. These further demonstrated crucial effect of miRNAs on skeletal muscle development.
In goats, the studies of miRNA expression profiles in skeletal muscle tissue have mainly been focused on different development stages. For example, Wang et al. identified 336 differentially expressed miRNAs in skeletal muscle of Huanghuai goats between fetal stage and 6-month-old stage, of which miR-424-5p, miR-29a, miR-129-3p, miR-181b, and miR-181d were involved in multiple important pathways related to muscle development (18). Guo et al. and Ling et al. also found some important differentially expressed miRNAs in skeletal muscle from prenatal stages to neonatal stage in Jianzhou Da′er goats and Anhui white goats (19, 20). However, little is known about the miRNA profiles of muscle tissues in other goat breeds, or between different goat breeds.
Ziwuling black (ZB) goats and Liaoning cashmere (LC) goats are both indigenous goat breeds in China. There are significant differences in meat production performance and muscle nutrients between the two breeds. For example, LC goats had higher carcass weight, muscle mass and intramuscular fat content, but lower muscle fiber density as well as contents of linoleic (C18: 2n-6), 11C, 14C-eicosadienoic acid (C20: 2n-6), moisture, and crude ash in meat when compared to ZB goats (P < 0.05) (21). In this context, elucidating the molecular mechanisms regulating these phenotypic differences between LC and ZB goats can provide insight for improving meat production performance of goats and other livestock. In this study, the expression profiles of miRNAs were compared in the Longissimus dorsi muscle between LC and ZB goats using small RNA sequencing. The differentially expressed miRNAs were screened between the two caprine breeds and the roles of miRNAs were also uncovered in skeletal muscle development and intramuscular fat deposition in goats.
Materials and Methods
Ethics Statement
All animal procedures in this study were approved by Animal Experiment Ethics Committee of Gansu Agricultural University with an approval number of GSAU-ETH-AST-2021-028.
Longissimus dorsi Muscle Sample Collection and RNA Extraction
Ten healthy, 9-month-old male goats (five LC goats and five ZB goats) were selected from Yongfeng Goat Breeding Company in Huan County, Gansu Province, China. All goats were raised under the same environmental conditions and nutrition levels. After being slaughtered, the Longissimus dorsi muscle samples from the area between 12th and 13th ribs on the left carcass of each goat were collected and then frozen in liquid nitrogen immediately until further use. The meat production performance, muscle fiber size and intramuscular fat content from these LC and ZB goats have been reported by Wang et al. (21) and were also presented in Supplementary Material 1.
Total RNA from Longissimus dorsi muscle samples was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). The concentration and purity of the RNA extracted were assessed using a Nanodrop 2000 (Thermo Scientific, MA, USA). Only samples with an RNA concentration >80 ng/uL and a purity of 1.80–2.10 were used for the study. The Agilent 2100 Bioanalyzer (Agilent, CA, USA) was used to assess RNA Integrity Number (RIN) of samples. Only RNA samples with RIN value ≥ 7 were used for small RNA enrichment.
Small RNA Library Construction and Sequencing
Ten small RNA libraries were generated using a TruSeq™ Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA) and then sequenced using an Illumina HiSeqTM4000 sequencer (Illumina, San Diego, CA, United States) at the Gene Denovo Biotechnology Co., Ltd (Guangzhou, China). The clean reads were obtained by removing the reads containing adapters, low quality reads with quality scores < Q20 (the proportion of read bases whose error rate is <1%) or with unknown nucleotides, and the reads shorter than 18nt in length in the raw reads, using fastp v0.18.0. First, the clean reads were mapped to GenBank database v209.0 and Rfam database v11.0 to annotate and remove other non-coding RNAs, including ribosomal RNA (rRNA), transfer RNA (tRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), and small cytoplasmic RNA (scRNA). Secondly, the clean reads were mapped to the Caprine Genome Assembly ARS1 (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/704/415/GCF_001704415.1_ARS1) to remove their exons, introns and repeated sequences. Subsequently, the remaining clean reads were searched against miRbase v22.0 to annotate known caprine miRNAs and known miRNAs from other species (named known species-conserved miRNA). Finally, for the reads that were not annotated to miRBase V22.0, but matched the Caprine Genome Assembly ARS1, they were used to predict novel miRNAs using the miReap v.0.2. To ensure the uniquely annotated results for the reads, the following annotation ranking was used: rRNA > caprine miRNA > caprine miRNA edit > species-conserved miRNA > repeat sequence > exon sequence > novel miRNA > intron sequence.
Differentially Expressed miRNAs Analysis and Small RNA Sequencing Results Validation
The expression level of miRNAs was first normalized using transcripts per million (TPM). The TPM value is calculated by actual reads of each miRNA*106 by total reads of all miRNAs. The DESeq v2.0 (22) was used to screen differentially expressed miRNAs in Longissimus dorsi muscle tissues between LC and ZB goats, using the thresholds of a |fold change| > 2.0 and P-value < 0.05. To validate the accuracy of small RNA sequencing results, 23 differentially expressed miRNAs were selected for reverse transcription-quantitative PCR (RT-qPCR) analysis, including eight up-regulated miRNAs (miR-628-5p, miR-885-3p, novel-m0312-3p, miR-1994-3p, miR-67-3p, miR-278-3p, miR-307-3p, and novel-m0298-5p) and 15 down-regulated miRNAs (miR-381, miR-127-3p, miR-200c, miR-136-3p, miR-487b-3p, miR-200a, miR-410-3p, miR-136-5p, miR-127-5p, miR-141, miR-200b, miR-276-3p, novel-m0213-5p, miR-2796-3p, and miR-429) in Longissimus dorsi muscle of LC goats compared to ZB goats. The same RNA samples as those used for the small RNA sequencing were used to generate cDNA using a miRNA 1st Strand cDNA Synthesis Kit (Accurate Biology, Hunan, China). The caprine U6 and 18SrRNA were used as internal references to normalize the relative expression level of miRNAs (18, 19). The RT-qPCR was performed in triplicate using 2 × ChamQ SYBR qPCR Master system (Vazyme, Nanjing, China) on an Applied Biosystems QuantStudio 6 Flex (Thermo Lifetech, MA, United States) platform. A 20 μL reaction system was used for the RT-qPCR analysis including 2.0 μL of the cDNA, 0.4 μL of each primer, 10 μL of SYBR qPCR master mix (Vazyme, Nanjing, China) and 7.2 μL of RNase-free water. The thermal profile included an initial denaturation of 30 s at 95°C, followed by 45 cycles of 95°C for 10 s, 60°C for 34 s, and 95°C for 15 s, and finished by 60°C for 60 s. The 2−ΔΔCt method was used to calculate the relative expression level of the miRNAs. The primer information used for RT-qPCR was presented in Supplementary Material 2.
Prediction, Validation, and Pathway Enrichment Analysis of the Target Genes of Differentially Expressed miRNAs
To investigate the potential roles of the differentially expressed miRNAs, miReap v0.2 (23), Miranda v3.3a (24), and TargetScan v7.0 (25) were used to predict their target genes and the predicted results from the three kinds of software were overlapped. To further verify the target relationship between the miRNAs and predicted target genes, a RT-qPCR analysis was performed to detect their relative expression levels in Longissimus dorsi muscle tissues of LC and ZB goats. The caprine GAPDH was used as an internal reference (5), and the primer information was listed in Supplementary Material 2. The RNA samples that were the same as those used for the small RNA sequencing analysis, were used to synthesize cDNA using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA, United States). The same conditions and thermal profiles described above were used to perform the RT-qPCR analysis. The Pearson's coefficients in expression levels between the miRNAs and the target genes were calculated using SPSS v24.0. For negatively correlative pairs of miRNA-mRNA, the Cytoscape v3.5.1 (26) was used to construct an interaction network. The enrichment analysis of the signaling pathway of the target genes was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (27). The significant pathways (P < 0.05) were defined by hypergeometric test, and the P-values were corrected using the calculated False Discovery Rate (FDR) value.
Results
Quality Control of Small RNA Sequencing Data
The concentration of 10 RNA samples collected from Longissimus dorsi muscle ranged from 124 to 272 ng/uL, while their purity ranged from 1.92 to 2.05 (Supplementary Material 3). On average, a total of 16,061,071 and 15,645,993 raw reads were generated from Longissimus dorsi muscle tissues of LC and ZB goats, respectively. The raw reads obtained in the study have been deposited in GenBank with accession numbers SRR16760528-SRR16760537. After removing low quality reads, adaptors and reads shorter than 18nt in length, an average of 15,337,813 and 14,815,683 clean reads were obtained in LC and ZB goats, respectively, of which 79.4 and 79.3% reads were mapped well to the caprine reference genome ARS1. Of these reads obtained, most of small RNA ranged from 18 to 24 nucleotides in length and the reads with 22 nucleotides were the most common, accounting for 45.7 and 45.4% of total reads in LC and ZB goats, respectively (Figure 1).
Figure 1. The nucleotide length distribution of small RNA reads obtained from Longissimus dorsi muscle tissues of Liaoning cashmere (LC) and Ziwuling black (ZB) goats.
Identification of Known and Novel Caprine miRNAs
A total of 1,623 miRNAs were detected in Longissimus dorsi muscle tissues from both LC and ZB goats, including 410 known caprine miRNAs, 928 known species-conserved miRNAs and 285 novel miRNAs (Supplementary Material 4). Of these miRNAs, 1,142 were co-expressed in both breeds, while 230 and 251 miRNAs were only expressed in LC and ZB goats, respectively. Among all the small RNAs annotated in this study, the known miRNAs including mature caprine miRNA and species-conserved miRNA were the most abundant, which represented 84.0 and 83.6% of the total number of small RNA reads in LC and ZB goats, respectively (Figure 2).
Figure 2. The percentage of small RNA types in Longissimus dorsi muscle of Liaoning cashmere (LC) and Ziwuling black (ZB) goats. Known miRNAs included mature caprine miRNAs and species-conserved miRNAs. Others RNAs included the sequences aligned to exon and intron, repeated sequences, caprine miRNA edit, and sequences that were not aligned to any database.
Of the 1,623 miRNAs identified, miR-133a-3p was the most abundant with an TPM value of 125,235 and 136,882 in LC goats and ZB goats, respectively, followed by miR-26a-5p, miR-1, miR-99a-5p, and miR-27b-3p (Supplementary Material 4). Most notably, miR-133a-3p and miR-1 are members of myomiRs (namely miRNAs specific to muscle tissues) (28). In addition, myomiRs also included other highly expressed miRNAs, such as miR-206, miR-133b, and miR-208b (28, 29).
Screening and Validation of Differentially Expressed miRNAs
A total of 159 miRNAs were identified to be differentially expressed in Longissimus dorsi muscle tissues when comparing LC goats and ZB goats. Twenty-four miRNAs had higher expression in LC goats compared to ZB goats including one known caprine miRNA, 16 known species-conserved miRNAs and seven novel miRNAs (Supplementary Material 5). Among these up-regulated miRNAs in LC goats, miR-1994-3p was the most significant differentially expressed miRNA, followed by miR-67-3p, miR-278-3p, miR-307-3p, and miR-503-5p.
A total of 135 down-regulated miRNAs were identified in LC goats including 40 known caprine miRNAs, 86 known species-conserved miRNAs and nine novel miRNAs (Supplementary Material 5). Of these miRNAs, the most prominent down-regulated miRNA was miR-381, followed by miR-276-3p, miR-429, miR-2796-3p, and miR-136-3p.
The results from the RT-qPCR for 23 differentially expressed miRNAs were in consistency with those obtained from the small RNA sequencing analysis (Figure 3). Because miR-1994-3p and miR-429 were only expressed in LC and ZB goats, respectively, the log2 fold-change for LC goats relative to ZB goats was infinity for these two miRNAs. In this context, their relative expression levels are not presented in Figure 3. These results demonstrate the repeatability and reliability of small RNA sequencing results.
Figure 3. RT-qPCR validation (A) of 21 differentially expressed miRNAs in Longissimus dorsi muscle tissues between Liaoning cashmere (LC) and Ziwuling black (ZB) goats identified using small RNA sequencing (B). The error bars represent standard deviation value of the means for three independent replicates for each sample.
Predication and KEGG Analysis of the Target Genes of Differentially Expressed miRNAs
The results from interaction analysis of miReap v0.2, Miranda v3.3a, and TargetScan v7.0 revealed a total of 15,029 target genes identified for the 159 differentially expressed miRNAs. To clearly exhibit the interaction between the miRNAs and their target genes, 12 differentially expressed miRNAs were further selected, including the four most up-regulated miRNAs (miR-1994-3p, miR-67-3p, miR-278-3p, and miR-307-3p) and the five most down-regulated miRNAs (miR-381, miR-276-3p, miR-429, miR-2796-3p, and miR-136-3p) in LC goats, as well as two novel up-regulated miRNAs (novel-m0312-3p and novel-m0298-5p) and one novel down-regulated miRNA (novel-m0213-5p) in LC goats. It was notable that the novel-m0312-3p and novel-m0298-5p had the highest expression levels in LC goats among all novel up-regulated miRNAs, while the novel-m0213-5p had the highest expression levels in ZB goats among all novel down-regulated miRNAs. There were 5,407 target genes in total for these 12 differentially expressed miRNAs, ranging from 87 target genes for miR-2796-3p to 1,295 targets for miR-429. For the 5,407 target genes, the genes related to muscle development and intramuscular fat deposition were further selected and their target relationships with corresponding miRNAs were verified using RT-qPCR (Figures 3, 4). As shown in Supplementary Material 6, there were negative correlations in expression levels between the 12 miRNAs selected and their target genes. These suggest potential target relationships between these miRNAs and their target genes. Finally, a miRNA-mRNA interaction network was constructed (Figure 5). Some functional genes that have been previously described to be related with skeletal muscle development and intramuscular fat deposition were identified in this analysis. For example, the target genes JAG2, IGFBP5, HDAC9, and FOXO1 were closely associated with myogenesis, while SOX6 and COL1A1 were reported to regulate adipogenesis (30–35) (Figure 5).
Figure 4. Relative expression level of the target genes of 12 miRNAs selected in Longissimus dorsi muscle between Liaoning cashmere (LC) and Ziwuling black (ZB) goats detected using RT-qPCR.
Figure 5. The miRNA-mRNA interaction network of 12 differentially expressed miRNAs and their target genes. The red triangles and inverted triangles represent up-regulated and down-regulated miRNAs in Longissimus dorsi muscle of Liaoning cashmere (LC) goats compared to Ziwuling black (ZB) goats, respectively. The green circles represent the target genes of the miRNAs.
To further investigate the possible function of the target genes of differentially expressed miRNAs identified, a KEGG pathway analysis was performed (Supplementary Material 7, Figure 6). For up-regulated miRNAs in Longissimus dorsi muscle of LC goats, the most enriched pathway with the lowest P-value was axon guidance (P-value = 2.55E-09), followed by Ras signaling pathway (P-value = 1.82E-07) and Rap1 signaling pathway (P-value = 4.57E-07; Figure 6A, Supplementary Material 7). In addition, FoxO signaling pathway was also found among the top 10 significant pathways (P-value = 4.46E-06; Figure 6A). For down-regulated miRNAs in LC goats, the target genes were most enriched in pathways in cancer (P-value = 3.77E-17), followed by Hippo signaling pathway (P-value = 1.17E-16) and metabolic pathways (P-value = 2.46E-16; Figure 6B, Supplementary Material 7).
Figure 6. The top 10 KEGG signaling pathways for the target genes of up-regulated (A) and down-regulated (B) miRNAs in the Longissimus dorsi muscle of Liaoning cashmere (LC) goats compared to Ziwuling black (ZB) goats. The left side Y-axis represents the number of the target genes of differentially expressed miRNAs involved in the pathway, while the Y axis on the right side shows the value of -Log10 (P-value).
Discussion
This study compared the Longissimus dorsi muscle tissue expression profiles of miRNAs of LC goats, with those of ZB goats that had lower carcass weight, muscle mass and intramuscular fat content. It was observed in the study that most miRNAs were 22 nucleotides in length identified in both caprine breeds and this was in accordance with the typical size range of mature miRNAs from Dicer-derived products (36). The results were also identical to the length distribution of small RNA reads in skeletal muscle tissues of other goat breeds (18, 20, 37), pigs (38), sheep (14, 15), and cattle (11). Additionally, our observation that the vast majority of small RNAs was known miRNAs in the two goat breeds has been also observed in previous studies of Longissimus dorsi muscle tissue (18, 20, 37). For example, Wang et al. found 81.5 and 82.5% of known miRNAs in muscle tissue of Huanghuai goats with fetal stage and 6-month-old stage, respectively (18). Similarly, 66.7 and 76.5% known miRNAs were found in skeletal muscle tissues of Anhui white goats (20) and Boer goats (37), respectively.
Of the top five highly expressed miRNAs found in both goat breeds, miR-133a-3p and miR-1 are members of myomiRs, which were specifically expressed in muscle tissues. The myomiRs have been reported to play key roles in regulating hypertrophy and regeneration of muscle fiber, as well as the differentiation and proliferation of muscle satellite cells (3, 39, 40). The miR-133a-3p was found to be the most highly expressed miRNA in Longissimus dorsi muscle of Boer goats (37). In Anhui white goats, Boer goats, Huanghuai goats and Jianzhou Da′er goats (18–20, 37), miR-1 was also one of the most abundant miRNAs in skeletal muscle tissues. The top five highly expressed miRNAs in the study also included miR-26a-5p, miR-99a-5p, and miR-27b-3p. The miR-99a-5p and miR-27b-3p have been reported to regulate proliferation and differentiation of skeletal muscle satellite cells in chicken (41) and goats (5), respectively, while miR-26a-5p promoted myogenesis by targeting Smad1 and Smad4 of TGF-β/BMP pathway (42). In short, these highly expressed miRNAs may be necessary for the growth and development of caprine skeletal muscle.
In this study, the expression levels of miR-200 family were very low only with 0-101.3 TPM values and also down-regulated in muscle tissues of LC goats, including miR-429, miR-200b, miR-200a, miR-200c, and miR-141. The miR-200 family acts primarily as a negative regulator of muscle cells and adipocytes development. For example, miR-200b suppressed proliferation of C2C12 myoblast (43) and differentiation of ovine preadipocytes (44). The miR-200c and miR-429 have been reported to inhibit differentiation of C2C12 myoblast (45), porcine preadipocytes, and C2C12 myoblast (46, 47), respectively. These suggest that down-regulated expression of miR-200 family may be responsible for the higher carcass weight, muscle fiber size and intramuscular fat content in LC goats.
The miR-381 was the most down-regulated miRNA in LC goats in the study. The differential expression of miR-381 in muscle tissues between different breeds with divergent meat production performance has also been found in sheep (15) and pigs (38). This may reflect the breed-specific expression pattern of miR-381. Although the molecular mechanism of miR-381 regulating on the growth and development of skeletal muscle is unclear, a miRNA-mRNA network showed that miR-381 would target some important functional genes, such as jagged canonical Notch ligand 2 (JAG2), insulin like growth factor binding protein 5 (IGFBP5), phosphatase and tensin homolog (PTEN), etc. (Figure 5). JAG2 and IGFBP5 are important components of notch signaling and IGF signaling pathways, respectively. The two signaling pathways promoted muscle growth and development by regulating the activity of muscle satellite cells (30, 31). It was inferred that higher expression level of miR-381 in ZB goats may result in lower carcass weight by more inhibiting the expression level of the target genes JAG2 and IGFBP5.
It was notable that known species-conserved miRNAs identified in this study may play key roles in the growth and development of caprine skeletal muscle, although their sequences have not been deposited in known caprine miRNA database. For example, up-regulated known species-conserved miR-885-3p in LC goats was reported to promote myoblasts proliferation in cattle (48). On the contrary, down-regulated species-conserved miR-370-3p in LC goats played a negative role in skeletal myogenesis in mice (49). The two species-conserved miRNAs may partly explain higher meat production performance in LC goats compared to ZB goats. As one of the most down-regulated miRNAs in LC goats, the known species-conserved miR-276-3p and miR-2796-3p would target SRY-box transcription factor 6 (SOX6; Figure 5), which promoted adipogenesis in human by activating adipogenic regulators including PPARγ, C/EBPα, and MEST (34). We therefore speculate that the down-regulation of the two species-conserved miRNAs in LC goats may promote adipogenesis by elevating expression of SOX6, resulting in increased deposition of intramuscular fat in LC goats. Similarly, the most up-regulated species-conserved miR-1994-3p may contribute to higher carcass weight and muscle mass of LC goats as it would target histone deacetylase 9 (HDAC9). The gene was found to play negative roles in muscle cell differentiation (32).
Four down-regulated known caprine miRNAs (miR-127-3p, miR-217-5p, miR-410-3p, and miR-487b-3p) in LC goats attracted our attention. The miR-127-3p has been reported to inhibit proliferation of C2C12 myoblast (50), as well as proliferation and differentiation of porcine skeletal muscle satellite cells (51). Similarly, miR-217-5p and miR-487b-3p inhibited the differentiation of skeletal muscle cells in rats (52) and mice (53), respectively. Additionally, the inhibition of miR-410-3p on the differentiation of adipocyte by targeting IRS-1 has also been reported in humans (54). It could be therefore inferred that the lower expression levels of the four miRNAs in LC goats may be responsible for its higher meat production performance and intramuscular fat content compared to ZB goats.
It was noteworthy that some novel miRNAs identified in this study were also differentially expressed between the two goat breeds. The roles of these miRNAs in the growth and development of muscle and adipose tissues may be reflected by the function of their target genes. For example, as a positive regulator of muscle atrophy (33), Forkhead box O1 (FOXO1) would be targeted by up-regulated novel-m0312-3p in LC goats. The up-regulation of novel-m0312-3p would result in higher muscle mass in LC goats by inhibiting the expression of FOXO1 and its effect on muscle atrophy. Meanwhile, down-regulated novel-m0213-5p found in LC goats would target SMAD family member 4 (SMAD4), which was contributed to the proliferation of porcine intramuscular preadipocyte (55). It was therefore inferred that down-regulation expression of novel-m0213-5p was responsible for higher intramuscular fat content in LC goats by less inhibition of SMAD4 in expression compared to ZB goats.
As might be expected, the target genes of differentially expressed miRNAs identified in the study were involved in the growth and development of skeletal muscle or adipose tissues. As one of the most enriched pathways for the target genes of up-regulated miRNAs in LC goats, Ras signaling has been reported to negatively regulate skeletal muscle myogenesis (56). Rap 1 signaling pathway negatively regulated adipocyte differentiation and was also associated with myogenic differentiation (57, 58). FoxO signaling pathway accelerated skeletal muscle atrophy by inducing proteolytic and apoptotic (59, 60). These enriched pathways found in the study have also been described previously. For example, Ras signaling pathway was significantly enriched by the target genes of differentially expressed miRNAs identified in skeletal muscle of yak with different ages (61), and Rap1 signaling pathway was enriched by differentially expressed genes in muscle tissues of pigs (62) and goats (63) during different development stages. The three pathways described above may partly explain why LC goats have higher carcass weight and intramuscular fat content compared to ZB goats. Hippo signaling pathway was one of the most enriched pathways for the target genes of down-regulated miRNAs in LC goats with larger muscle fiber size and carcass weight. This is not surprising as the pathway is necessary for the increase of skeletal muscle mass (64, 65). Interestingly, MAPK signaling pathway was significantly enriched by both the target genes of up-regulated miRNAs and down-regulated miRNAs in LC goats. This suggests that the pathway may play dual roles in muscle development. The speculation was subsequently confirmed. MAPK signaling pathway has been mainly recognized as a positive regulator of myogenesis in animals (66). However, an inhibition effect on myogenesis in mice was also reported by Weston et al. (67).
Conclusion
This study compared the skeletal muscle tissue expression profiles of miRNAs between different goat breeds. The miR-381, miR-127-3p, miR-200b, miR-200c, miR-429, miR-217-5p, miR-885-3p, miR-370-3p, miR-1994-3p, miR-487b-3p, and novel-m0312-3p were found to be associated with muscle development, while miR-200b, miR-429, miR-276-3p, miR-2796-3p, miR-410-3p, and novel-m0213-5p were related to intramuscular fat deposition in goats.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in GenBank [accession: SRR16760528-SRR16760537].
Ethics Statement
The animal study was reviewed and approved by Animal Experiment Ethics Committee of Gansu Agricultural University (Ethic approval file No. GSAU-ETH-AST-2021-028). Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
JS and JW did the data analysis and wrote the manuscript. ZH, YuzL, HZ, YaL, JH, and XL performed investigation and collected the samples. SL, ZZ, YuaL, SY, and LW performed the formal analysis, methodology, and software. JW did the project administration and revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was funded by the fund for Basic Research Creative Groups of Gansu Province (18JR3RA190), the Fuxi Young Talents Fund of Gansu Agricultural University (Gaufx-02Y02), and the Projects of Gansu Agricultural University (GSAU-ZL-2015-033).
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/fvets.2022.911166/full#supplementary-material
References
1. Carrington JC, Ambros V. Role of microRNAs in plant and animal development. Science. (2003) 301:336–8. doi: 10.1126/science.1085242
3. Chen JF, Tao Y, Li J, Deng Z, Yan Z, Xiao X, et al. MicroRNA-1 and microRNA-206 regulate skeletal muscle satellite cell proliferation and differentiation by repressing Pax7. J Cell Biol. (2010) 190:867–79. doi: 10.1083/jcb.200911036
4. Hitachi K, Nakatani M, Tsuchida K. Myostatin signaling regulates Akt activity via the regulation of miR-486 expression. Int J Biochem Cell B. (2014) 47:93–103. doi: 10.1016/j.biocel.2013.12.003
5. Ling YH, Sui MH, Zheng Q, Wang KY, Wu H, Li WY, et al. MiR-27b regulates myogenic proliferation and differentiation by targeting Pax3 in goat. Sci Rep. (2018) 8:3909. doi: 10.1038/s41598-018-22262-4
6. Lu J, Ye H, Hui W, Miao Y, Li X, Cao J, et al. Transcriptome analysis of mRNA and miRNA in skeletal muscle indicates an important network for differential Residual Feed Intake in pigs. Sci Rep. (2015) 5:11953. doi: 10.1038/srep11953
7. Sun J, Xie M, Huang Z, Li H, Chen T, Sun R, et al. Integrated analysis of non-coding RNA and mRNA expression profiles of 2 pig breeds differing in muscle traits. J Anim Sci. (2017) 95:1092. doi: 10.2527/jas2016.0867
8. Wang Q, Qi R, Wang J, Huang W, Wu Y, Huang X, et al. Differential expression profile of miRNAs in porcine muscle and adipose tissue during development. Gene. (2017) 618:49–56. doi: 10.3390/genes8100271
9. Wei W, Li B, Liu K, Jiang A, Dong C, Jia C, et al. Identification of key microRNAs affecting drip loss in porcine longissimus dorsi by RNA-Seq. Gene. (2018) 647:276–82. doi: 10.1016/j.gene.2018.01.005
10. Huang YZ, Sun JJ, Zhang LZ, Li CJ, Womack JE, Li ZJ, et al. Genome-wide DNA methylation profiles and their relationships with mRNA and the microRNA transcriptome in bovine muscle tissue (Bos taurine). Sci Rep. (2014) 4:6546. doi: 10.1038/srep06546
11. Sun J, Sonstegard TS, Li C, Huang Y, Li Z, Lan X, et al. Altered microRNA expression in bovine skeletal muscle with age. Anim Genet. (2015) 46:227–38. doi: 10.1111/age.12272
12. Li N, Zhang Y, Li HP, Han L, Yan XM, Li HB, et al. Differential expression of mRNA-miRNAs related to intramuscular fat content in the longissimus dorsi in Xinjiang brown cattle. PLoS ONE. (2018) 13:e0206757. doi: 10.1371/journal.pone.0206757
13. Kappeler B, Regitano L, Poleti M, Cesar A, Moreira G, Gasparin G, et al. MiRNAs differentially expressed in skeletal muscle of animals with divergent estimated breeding values for beef tenderness. BMC Mol Biol. (2019) 20:1. doi: 10.1186/s12867-018-0118-3
14. Liu Z, Li C, Li X, Yao Y, Ni W, Zhang X, et al. Expression profiles of microRNAs in skeletal muscle of sheep by deep sequencing. Asian Austral J Anim. (2019) 32:757–66. doi: 10.5713/ajas.18.0473
15. Sun L, Lu S, Bai M, Xiang L, Li J, Jia C, et al. 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 Cell Biol. (2019) 38:423–35. doi: 10.1089/dna.2018.4408
16. Kaur M, Kumar A, Siddaraju NK, Fairoze MN, Chhabra P, Ahlawat S, et al. Differential expression of miRNAs in skeletal muscles of Indian sheep with diverse carcass and muscle traits. Sci Rep. (2020) 10:16332. doi: 10.1038/s41598-020-73071-7
17. Yuan C, Zhang K, Yue Y, Guo T, Liu J, Niu C, et al. Analysis of dynamic and widespread lncRNA and miRNA expression in fetal sheep skeletal muscle. Peer J. (2020) 8:e9957. doi: 10.7717/peerj.9957
18. Wang Y, Zhang C, Fang X, Zhao Y, Chen X, Sun J, et al. Identification and profiling of microRNAs and their target genes from developing caprine skeletal muscle. PLoS ONE. (2014) 9:e96857. doi: 10.1371/journal.pone.0096857
19. Guo J, Wei Z, Zhan S, Li L, Zhong T, Wang L, et al. Identification and expression profiling of miRNAome in goat longissimus dorsi muscle from prenatal stages to a neonatal stage. PLoS ONE. (2016) 11:e0165764. doi: 10.1371/journal.pone.0165764
20. Ling Y, Zheng Q, Jing J, Sui M, Zhu L, Li Y, et al. RNA-Seq reveals miRNA role shifts in seven stages of skeletal muscles in goat fetuses and kids. Front Genet. (2020) 11:684. doi: 10.3389/fgene.2020.00684
21. Wang JQ, Shen JY, Liu X, Li SB, Luo YZ, Zhao ML, et al. Comparative analysis of meat production traits, meat quality, and muscle nutrient and fatty acid contents between Ziwuling black goats and Liaoning cashmere goats. Acta Pratacult Sin. (2021) 30:166–77. doi: 10.11686/cyxb2020199
22. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
23. Liu W, Xu L, Wang Y, Shen H, Zhu X, Zhang K, et al. Transcriptome-wide analysis of chromium-stress responsive microRNAs to explore miRNA-mediated regulatory networks in radish (Raphanus sativus L). Sci Rep. (2015) 5:14024. doi: 10.1038/srep14024
24. Turner DA. Miranda: a non-strict functional language with polymorphic types. In: Conference on Functional Programming Languages and Computer Architecture. Berlin; Heidelberg: Springer (1985). p. 1–16. doi: 10.1007/3-540-15975-4_26
25. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. (2005) 120:15–20. doi: 10.1016/j.cell.2004.12.035
26. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
27. Minoru K, Michihiro A, Susumu G, Masahiro H, Mika H, Masumi I, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. (2008) 36:480–4. doi: 10.1093/nar/gkm882
28. McCarthy JJ. The MyomiR network in skeletal muscle plasticity. Exerc Sport Sci Rev. (2011) 39:150–4. doi: 10.1097/JES.0b013e31821c01e1
29. Li D, Xia L, Chen M, Lin C, Wu H, Zhang Y et al. miR-133b, a particular member of myomiRs, coming into playing its unique pathological role in human cancer. Oncotarget. (2017) 8:50193–208. doi: 10.18632/oncotarget.16745
30. Shan T, Xu Z, Wu W, Liu J, Wang Y. Roles of Notch1 signaling in regulating satellite cell fates choices and postnatal skeletal myogenesis. J Cell Physiol. (2017) 232:2964–7. doi: 10.1002/jcp.25730
31. Zhang L, Wang XH, Wang H, Du J, Mitch WE. Satellite cell dysfunction and impaired IGF-1 signaling cause CKD-induced muscle atrophy. J Am Soc Nephrol. (2010) 21:419–27. doi: 10.1681/ASN.2009060571
32. Haberland M, Arnold MA, McAnally J, Phan D, Kim Y, Olson EN. Regulation of HDAC9 gene expression by MEF2 establishes a negative-feedback loop in the transcriptional circuitry of muscle differentiation. Mol Cell Biol. (2007) 27:518–25. doi: 10.1128/MCB.01415-06
33. Oyabu M, Takigawa K, Mizutani S, Hatazawa Y, Fujita M, Ohira Y, et al. FOXO1 cooperates with C/EBPδ and ATF4 to regulate skeletal muscle atrophy transcriptional program during fasting. FASEB J. (2022) 36:e22152. doi: 10.1096/fj.202101385RR
34. Leow SC, Poschmann J, Too PG, Yin J, Joseph R, McFarlane C, et al. The transcription factor SOX6 contributes to the developmental origins of obesity by promoting adipogenesis. Development. (2016) 143:950–61. doi: 10.1242/dev.131573
35. Côté JA, Guénard F, Lessard J, Lapointe M, Biron S, Vohl MC, et al. Temporal changes in gene expression profile during mature adipocyte dedifferentiation. Int J Genomics. (2017) 2017:5149362. doi: 10.1155/2017/5149362
36. Zhang B, Stellwag EJ, Pan X. Large-scale genome analysis reveals unique features of microRNAs. Gene. (2009) 443:100–9. doi: 10.1016/j.gene.2009.04.027
37. Ling YH, Ding JP, Zhang XD, Wang LJ, Zhang YH, Li YS, et al. Characterization of microRNAs from goat (Capra hircus) by Solexa deep-sequencing technology. Genet Mol Res. (2013) 12:1951–61. doi: 10.4238/2013.June.13.4
38. Hou X, Yang Y, Zhu S, Hua C, Zhou R, Mu Y, et al. Comparison of skeletal muscle miRNA and mRNA profiles among three pig breeds. Mol Genet Genomics. (2016) 291:559–73. doi: 10.1007/s00438-015-1126-3
39. Chen JF, Mandel EM, Thomson JM, Wu Q, Callis TE, Hammond SM, et al. The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation. Nat Genet. (2006) 38:228–33. doi: 10.1038/ng1725
40. Nakasa T, Ishikawa M, Shi M, Shibuya H, Adachi N, Ochi M. Acceleration of muscle regeneration by local injection of muscle-specific microRNAs in rat skeletal muscle injury model. J Cell Mol Med. (2010) 14:2495–505. doi: 10.1111/j.1582-4934.2009.00898.x
41. Cao X, Tang S, Du F, Li H, Shen X, Li D, et al. MiR-99a-5p regulates the proliferation and differentiation of skeletal muscle satellite cells by targeting MTMR3 in chicken. Genes. (2020) 11:369. doi: 10.3390/genes11040369
42. Dey BK, Gagan J, Yan Z, Dutta A. MiR-26a is required for skeletal muscle differentiation and regeneration in mice. Gene Dev. (2012) 26:2180–91. doi: 10.1101/gad.198085.112
43. Yao CX, Wei QX, Zhang YY, Wang WP, Xue LX, Yang F, et al. MiR-200b targets GATA-4 during cell growth and differentiation. RNA Biol. (2013) 10:465–80. doi: 10.4161/rna.24370
44. Jin X, Wang J, Hu J, Liu X, Li S, Lu Y, et al. MicroRNA-200b regulates the proliferation and differentiation of ovine preadipocytes by targeting p27 and KLF9. Animals. (2021) 11:2417. doi: 10.3390/ani11082417
45. D'Agostino M, Torcinaro A, Madaro L, Marchetti L, Sileno S, Beji S, et al. Role of miR-200c in myogenic differentiation impairment via p66shc: implication in skeletal muscle regeneration of dystrophic mdx mice. Oxid Med Cell Longev. (2018) 2018:4814696. doi: 10.1155/2018/4814696
46. Peng Y, Chen FF, Ge J, Zhu JY, Shi XE, Li X, et al. MiR-429 inhibits differentiation and promotes proliferation in porcine preadipocytes. Int J Mol Sci. (2016) 17:2047. doi: 10.3390/ijms17122047
47. Nguyen MT, Min KH, Lee W. Palmitic acid-induced miR-429-3p impairs myoblast differentiation by downregulating CFL2. Int J Mol Sci. (2021) 22:10972. doi: 10.3390/ijms222010972
48. Elsaeid Elnour I, Dong D, Wang X, Zhansaya T, Khan R, Jian W, et al. Bta-miR-885 promotes proliferation and inhibits differentiation of myoblasts by targeting MyoD1. J Cell Physiol. (2020) 235:6625–36. doi: 10.1002/jcp.29559
49. Zhang P, Du J, Guo X, Wu S, He J, Li X, et al. LncMyoD promotes skeletal myogenesis and regulates skeletal muscle fiber-type composition by sponging miR-370-3p. Genes. (2021) 12:589. doi: 10.3390/genes12040589
50. Yuan R, Zhang X, Fang Y, Nie Y, Cai S, et al. MiR-127-3p inhibits the proliferation of myocytes by targeting KMT5a. Biochem Bioph Res Co. (2018) 503:970–6. doi: 10.1016/j.bbrc.2018.06.104
51. Chen X, Zhao C, Dou M, Sun Y, Yu T, et al. Deciphering the miRNA transcriptome of Rongchang pig longissimus dorsi at weaning and slaughter time points. J Anim Physiol An N. (2020) 104:954–64. doi: 10.1111/jpn.13314
52. Zhu M, Chen G, Yang Y, Yang J, Qin B, Gu L. MiR-217-5p regulates myogenesis in skeletal muscle stem cells by targeting FGFR2. Mol Med Rep. (2020) 22:850–8. doi: 10.3892/mmr.2020.11133
53. Wang J, Tan J, Qi Q, Yang L, Wang Y, Zhan C, et al. MiR-487b-3p suppresses the proliferation and differentiation of myoblasts by targeting IRS1 in skeletal muscle myogenesis. Int J Biol Sci. (2018) 14:760–74. doi: 10.7150/ijbs.25052
54. Sun D, Ding Z, Shen L, Yang F, Han J, Wu G. MiR-410-3P inhibits adipocyte differentiation by targeting IRS-1 in cancer-associated cachexia patients. Lipids Health Dis. (2021) 20:115. doi: 10.1186/s12944-021-01530-9
55. Zhang Q, Cai R, Tang G, Zhang W, Pang W. MiR-146a-5p targeting SMAD4 and TRAF6 inhibits adipogenensis through TGF-β and AKT/mTORC1 signal pathways in porcine intramuscular preadipocytes. J Anim Sci Biotechno. (2021) 12:12. doi: 10.1186/s40104-020-00525-3
56. Mitin N, Kudla AJ, Konieczny SF, Taparowsky EJ. Differential effects of Ras signaling through NFkappaB on skeletal myogenesis. Oncogene. (2001) 20:1276–86. doi: 10.1038/sj.onc.1204223
57. Pizon V, Cifuentes-Diaz C, Mège R, Baldacci G, Rieger F. Expression and localization of RAP1 proteins during myogenic differentiation. Eur J Cell Biol. (1996) 69:224.
58. Yeung F, Ramírez CM, Mateos-Gomez PA, Pinzaru A, Ceccarini G, Kabir S, et al. Nontelomeric role for Rap1 in regulating metabolism and protecting against obesity. Cell Rep. (2013) 3:1847–56. doi: 10.1016/j.celrep.2013.05.032
59. Mcloughlin TJ, Smith SM, Delong AD, Wang H, Unterman TG, Esser KA. FoxO1 induces apoptosis in skeletal myotubes in a DNA-binding-dependent manner. Am J Physiol Cell Physiol. (2009) 297:C548–555. doi: 10.1152/ajpcell.00502.2008
60. Egerman MA, Glass DJ. Signaling pathways controlling skeletal muscle mass. Crit Rev Biochem Mol Biol. (2014) 49:59–68. doi: 10.3109/10409238.2013.857291
61. Ji H, Wang H, Ji Q, Ji W, Zhong J. Differential expression profile of microRNA in yak skeletal muscle and adipose tissue during development. Genes Genom. (2020) 42:1347–59. doi: 10.1007/s13258-020-00988-8
62. Wang Y, Wang J, Hu H, Wang H, Wang C, Lin H, et al. Dynamic transcriptome profiles of postnatal porcine skeletal muscle growth and development. BMC Genom Data. (2021) 22:32. doi: 10.1186/s12863-021-00984-1
63. Zhan S, Zhao W, Song T, Dong Y, Zhang H. Dynamic transcriptomic analysis in hircine longissimus dorsi muscle from fetal to neonatal development stages. Funct Integr Genomic. (2018) 18:43–54. doi: 10.1007/s10142-017-0573-9
64. Watt KI, Turner BJ, Hagg A, Zhang X, Davey JR, Qian H, et al. The Hippo pathway effector YAP is a critical regulator of skeletal muscle fiber size. Nat Commun. (2015) 6:6048. doi: 10.1038/ncomms7048
65. Olouyomi G, Marc F, Louise D. Hippo pathway and skeletal muscle mass regulation in mammals: a controversial relationship. Front Physiol. (2017) 8:190. doi: 10.3389/fphys.2017.00190
66. Keren A, Tamir Y, Bengal E. The p38 MAPK signaling pathway: a major regulator of skeletal muscle development. Mol Cell Endocrinol. (2006) 252:224–30. doi: 10.1016/j.mce.2006.03.017
Keywords: microRNA (miRNA), muscle development, intramuscular fat, small RNA sequencing, goat
Citation: Shen J, Hao Z, Luo Y, Zhen H, Liu Y, Wang J, Hu J, Liu X, Li S, Zhao Z, Liu Y, Yang S and Wang L (2022) Deep Small RNA Sequencing Reveals Important miRNAs Related to Muscle Development and Intramuscular Fat Deposition in Longissimus dorsi Muscle From Different Goat Breeds. Front. Vet. Sci. 9:911166. doi: 10.3389/fvets.2022.911166
Received: 02 April 2022; Accepted: 09 May 2022;
Published: 13 June 2022.
Edited by:
Li Ma, University of Maryland, United StatesReviewed by:
Zeying Wang, Shenyang Agricultural University, ChinaRui Su, Inner Mongolia Agricultural University, China
Copyright © 2022 Shen, Hao, Luo, Zhen, Liu, Wang, Hu, Liu, Li, Zhao, Liu, Yang and Wang. 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.
*Correspondence: Jiqing Wang, wangjq@gsau.edu.cn