RNA-Seq Reveals miRNA Role Shifts in Seven Stages of Skeletal Muscles in Goat Fetuses and Kids

MicroRNAs (miRNAs) are indispensable for the regulation of skeletal muscle. We performed RNA sequencing (RNA-seq) to establish a comprehensive miRNA profiling of goats in seven stages, namely, 45- (F45), 65- (F65), 90- (F90), 120- (F120), and 135-day (F135) fetuses, newborn (B1), and 90-day-old (B90) kids. In total, 421 known miRNAs and 228 goat novel miRNAs were identified in the data, and the average abundance of 19 miRNAs in seven stages exceeds 10,000 reads per million. Furthermore, 420 differentially expressed miRNAs (DEmiRNAs) were identified in all comparison group at seven stages, 80 of which were uniquely differentially expressed in the B1 and B90 comparison groups. Pathway analysis indicated that this group was associated with the release of muscle hypertrophy and regulation of myoblast proliferation. Besides, 305 DEmiRNAs were clustered into three significantly enriched profiles (profiles 11, 16, and 19). Function analysis revealed that profile 16 was related to muscle hypertrophy and differentiation. Profile 11 was involved in multiple enzyme activities and metabolic processes in muscle cells. And profile 19 was involved in material transport and structural stability. Two highly expressed miRNAs and three key miRNAs (chi-miR-328-3p, chi-miR-767, and chi-miR-150) of these profiles were verified to be consistent with the data by quantitative real-time PCR. These results provided a catalog of goat muscle-associated miRNAs, allowing us to better understand the transformation of miRNA roles during mammalian muscle development.


INTRODUCTION
Skeletal muscle is a complex and important tissue that occupies 40% of mammals' body weight (Jung et al., 2019). The development of embryonic skeletal muscle tissue is mainly manifested by an increase in the number of muscle fibers. The myoblasts from the mesoderm first migrate to the limbs and trunk and rapidly proliferate (Endo, 2015;Chal and Pourquie, 2017). After reaching a certain number, myoblasts differentiate into myocytes and fuse into multinucleated myotubes (Endo, 2015;Ganassi et al., 2018). The fetal period is the main period of muscle formation, and the number of muscle fibers remains constant after birth (Zhu et al., 2004;Wang et al., 2010). After the mammal is born, the growth of skeletal muscle tissue mainly depends on the circumference and length of the muscle fiber, also known as muscle fiber hypertrophy (Al-Shanti and Stewart, 2009;Murach et al., 2018). Therefore, there are significant differences in the development of skeletal muscles before and after birth. And the transformation of skeletal muscle growth and development is essentially the result of gene regulation. MicroRNAs (miRNAs), the non-coding single-stranded RNA of 19-25 nucleotides, are emerging as robust players of gene regulation (Adlakha and Saini, 2014). Most of the miRNAs in animals mainly bind to the 3'-UTR, 5'-UTR, or CDS sequence of their target mRNA through their 5'-end twoto eight-seed sequences, resulting in interference or degradation of the target mRNA (Meister, 2007;Helwak et al., 2013). Thus, they are involved in various biological processes, such as cell proliferation (Roy et al., 2017), apoptosis (Gorbea et al., 2017), differentiation (Title et al., 2018), metabolism (Antoniali et al., 2017), and tumor metastasis (Fang et al., 2018). And miRNAs have different expression patterns before and after animal birth in skeletal muscle . Analysis of this difference will help in the further study of the regulation mechanism of miRNAs on animal skeletal muscle development.
Goat, one of the most basic commercial animals, has attracted increasing attention as a viable source of meat. However, goat has lower meat production than other commercial animals, which hinders the development of goat industries. Some studies have shown that miRNAs have an indispensable role in goat muscle development at different periods. For example, miR-27b was identified to inhibit the proliferation but promoted differentiation of skeletal muscle satellite via targeting PAX3 (Ling et al., 2018). And miR-487b-3p, which is significantly higher in the fetus than in adult, suppressed myoblast proliferation and differentiation by targeting IRS1 in skeletal muscle myogenesis .
In our study, seven developmental stages of Anhui white goat skeletal muscle were used for RNA sequencing (RNA-seq) analysis, namely, 45-(F45), 65-(F65), 90-(F90), 120 (F120), and 135-day (F135) fetuses and kids born within 24 h and 90 days (B1 and B90). miRNAs (known miRNAs and goat novel miRNAs) were systematically identified from multiple types of small RNAs (sRNA), and functions and pathways of dynamic expression profiles were elucidated. These findings revealed the molecular roles of miRNAs and the mechanisms operating at different stages of mammalian skeletal muscle development.

Sample Collection
A group of female Anhui white goats managed and raised in the Hefei Boda Animal Husbandry Technology Development Co., Ltd (Hefei, Anhui, China), was selected and treated with EAZI-Breed CIDR (CIDR, Hamilton, New Zealand) for 12 days for simultaneous estrus. The female goats with successful estrus were selected for artificial insemination. Fertilization was observed regularly. Twenty-one pregnant goats were used to obtain samples, each stage with three goats. Before sample collection, all the experimental animals (including female goats, fetuses, and kids) were injected with Jingsongling (2,4-xylyl xylazole, lot number 030725, produced by Shandong Zibo Veterinary Medicine Factory, Shandong, China) in the hips at a dose of 2 mg kg −1 for deep anesthesia. Under complete anesthesia, all experimental samples were killed by arterial bleeding with the help of slaughterhouse professionals (this method complied with the euthanasia guide of the Chinese Society of Laboratory Animals, no. T/CALAS 31-2017). Fetuses were obtained from goats that had been pregnant for 45, 65, 90, 120, and 135 days by cesarean section. Kids that are 1-and 90-day-old were also collected. The longissimus dorsi muscles of the seven stages of goat development were collected as experimental samples. Each stage had three biological replicates. After collection, the samples were washed three times with PBS and placed in liquid nitrogen to be stored frozen.

Library Preparation and Sequencing of Skeletal Muscle
Total RNA isolation was performed using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) on 21 longissimus dorsi muscle samples. Then, RNA concentration was measured using the Qubit R RNA Assay Kit in Qubit R 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States), RNA degradation and contamination were checked with 1% agarose gels, and RNA integrity was assessed using the RNA Nano 6000 Assay Kit (Agilent Technologies, Santa Clara, CA, United States) of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). A total amount of 3 µg total RNA per sample was used as input material for the sRNA sequencing libraries generated using NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB, Ipswich, MA, United States), and index codes were added to attribute sequences to each sample, followed by reverse transcription. PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for Illumina, and index (X) primer. PCR products were purified on an 8% polyacrylamide gel (100 V, 80 min). The small non-coding RNA plus the 3' and 5' adaptors were 140-160 bp in length, and the gel of these fractions was recovered and dissolved in 8 µl of elution buffer. Then, the quality of libraries was calculated by a DNA high-sensitivity chip on the Agilent Bioanalyzer 2100 system. Finally, the clustering of the indexcoded samples was performed on a cBot Cluster Generation System using a TruSeq SR Cluster Kit v3-cBot-HS (Illumina), and the library preparations were sequenced on an Illumina HiSeq 2500 platform, and 50-bp single-end reads were generated.

Data Calculation of Small RNA
Raw data in a fastq format were first processed through custom Perl and Python scripts. Raw data filtered reads contain ploy-N (with 5' adapter contaminants and ploy A or T or G or C, without 3' adapter or the insert tag) and low-quality reads to obtain clean data. The mismatched miRNA tags were mapped to the reference sequence 1 by Bowtie (Langmead et al., 2009). Then, the 18-to 35-nt fragments were selected from clean data for subsequent analyses. And the expression and distribution were analyzed according to the reference. Some sRNA reads map to multiple categories. In order to map each sRNA to only one annotation, we set the following precedence rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > novel miRNA.

Known and Goat Novel miRNA Alignments and Annotations
Small RNA reads were mapped to obtain known miRNAs. The total rRNA proportion was used as a marker as sample quality indicator. The acquisition of potential miRNAs and the mapping of secondary structures were performed by the modification software miRDeep2 (Friedlander et al., 2012) with miRBase 21.0 as a reference. sRNA tags were mapped to RepeatMasker, and the Rfam database was used to remove tags originating from protein-coding genes (intro ± link and exon ± link), repeat sequences, rRNA, tRNA, snRNA, and snoRNA. Moreover, the available software miREvo (Wen et al., 2012) and miRDeep2 (Friedlander et al., 2012) explored the secondary structure of unannotated fragments. And the Dicer cleavage site also has the minimum free energy of sRNA tags to predict goat novel miRNAs in the previous steps.
Analysis of miRNA miRNA families were identified from samples in other species. In our analysis pipeline, goat novel miRNA precursors were submitted to Rfam 2 to look for Rfam families.
Reads per million (RPM) was used as a unit for estimating the level of miRNA (Zhou et al., 2010). Differential expression analysis of two stages was performed by the DESeq R package (1.8.3). The P-value was adjusted using the Benjamini-and-Hochberg method. A corrected P-value (P-adj) < 0.05 of miRNA was considered a differential expression. Bubble charts were constructed using the OmicShare platform for data analysis 3 . And heat map analysis was formed using Morpheus 4 .

Functional and Pathway Analysis of miRNA Targets
The target genes of miRNAs were predicted using miRanda. Function and pathway differentially expressed miRNA (DEmiRNAs) target genes were demonstrated using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). GOseq-based Wallenius non-central hypergeometric distribution, which could be adjusted for gene length bias, was implemented for GO analysis (Young et al., 2010). KOBAS software was used to calculate candidate enrichment of target genes in the KEGG pathway (Mao et al., 2005).

Quantitative Real-Time PCR
Reverse transcription of selected DEmiRNAs was performed using a GoScript Reverse Transcriptase kit (Promega, Madison, WI, United States) according to the manufacturer's guide. All of the primer pairs (Table 1) were designed using NCBI and synthesized by the Shanghai General Biotech Co. Ltd (Table 1). U6 was amplified as a control. Next, GoTaq qPCR Master Mix (Promega, catalog number: A6002, Madison, WI, United States) was used to perform quantitative real-time PCR (qRT-PCR) on a LightCycler 96 (ABI) real-time PCR instrument according to the manufacturer's instructions. The PCR mixture included 7.5 µl SYBR Real-Time PCR System, 1.5 µl cDNA, and 1 µl forward and reverse primers (total 10 µM). The reaction was denatured at 95 • C for 2 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 1 min, and an extension at 60 • C for 30 s, and then recovered at a rate of 0.2 • C per 1 s to 95 • C; finally, it was stored at room temperature for 10 s. Each sample was run in triplicate. The target sequence quantity was normalized to the reference sequence and calculated as 2 − Ct . The Windows version of SPSS 21.0 5 was used for statistical analysis of standardized data. Data were presented as means ± SE and were considered statistically significant if the P-value was <0.05.

Overview of Seven-Stage Skeletal Muscle miRNA Sequencing Data
To identify miRNAs involved in early skeletal muscle development, sRNA libraries were constructed from total RNA of 21 individuals in seven stages. The libraries generated in seven developmental stages ranged from 10,217,979 to 15,858,768 raw reads, with an average of 13,409,245 reads per stage. After discarding adaptors and low-quality reads, an average of 13,180,785 clean reads was obtained per sample of the 21 sRNA libraries (Supplementary Table S1). After length 5 https://www.ibm.com/support/pages/spss-statistics-210-available-download  Table S1). Subsequently, the annotated sRNAs varying from 18 to 35 nt were divided into several categories: known miRNA (66.6%), intro (+ and -link, 2.04%), exon (+ and -link, 1.28%), rRNA (0.68%), snoRNA (0.20%), novel miRNA (0.25%), snRNA (0.02%), and others (28.19%) ( Figure 1A). Among them, known miRNAs accounted for the largest proportion. In addition, the obtained sRNAs were counted in length. Their lengths were concentrated at 20-24 nt ( Figure 1B). These data validated the reliability of the experiment. And there was a high correlation within samples at each stage, which also verified the reliability between samples at the same stage ( Figure 1C). A total of 421 known miRNAs and 228 goat novel miRNAs were obtained in the data, and their lengths were concentrated between 21 and 23 nt (Supplementary Table S2 and Figure 2A). Due to the specificity of the cleavage site, the first base of the miRNA mature sequence had a strong bias. And the main base preferences were U and A, which also verified the reliability of the data ( Figure 2B). Principal component analysis (PCA) of all detected miRNA showed that the developmental order could be accurately captured from F45 to B90, and there was a segmentation before and after birth ( Figure 2C). Pearson correlation analysis of all 21 sample pairs demonstrated similar results ( Figure 1C). These two calculations showed the same cluster expression profiles for the identified miRNAs.

Transcriptional Dynamics of Differentially Expressed miRNA
To understand the functional transformation of miRNAs at different stages in the skeletal muscle, it is essential to study the expression profiles of DEmiRNAs in skeletal muscle development at different stages. A total of 420 DEmiRNAs (P-adj < 0.05, including 364 known and 56 novel DEmiRNAs) were obtained in comparison groups of any two phases (Supplementary Table S3). The unbiased hierarchical clustering of DEmiRNAs indicated that the largest segmentation occurred between F90 and F120 ( Figure 3A). Moreover, in the comparison groups of two consecutive stages, F90 and F120 comparison groups also have the most DEmiRNAs ( Figure 3B). There were 76 DEmiRNAs observed between F90 and F120, of which 38 were upregulated and 38 were downregulated ( Figure 3B). Besides, the maximum number of DEmiRNAs (300) found between F45 and B90 and the minimum number of DEmiRNAs (5) were observed between F65 and F90, including chi-miR-136-5p, chi-miR-543-3p, chi-miR-497-5p, chi-miR-15b-3p, and chi-miR-34b-3p ( Figure 3B).

Stage-Specific DEmiRNA of Developmental Transition
To determine the role of DEmiRNAs throughout the seven developmental stages, we identified DEmiRNAs that were highly expressed at all stages. The total expression level of DEmiRNAs in seven stages varied between 0.14 and 1,632,833.52 with an intermediate value of 171.26 (Supplementary Table S3). The classification of small RNAs for goat skeletal muscle. Small RNA classification was performed using the priority of "known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > gene > novel miRNA." Others represent samples aligned to the reference sequence, but no alignment to the small RNA shown in the figure. (B) Percentage distribution of all small RNA sequence lengths. Characteristics of known and goat novel miRNAs in goat skeletal muscle. (C) Heat map of replicate samples of the same stage of goat skeletal muscle. The color spectrum, ranging from white to blue, represents Pearson correlation coefficients ranging from 1 to 0.7, indicating high to low correlations.
Frontiers in Genetics | www.frontiersin.org  Among them, 33.73% (142) DEmiRNAs showed high-abundance expression in seven stages (mean total RPM > 1,000), and the average expression level of 19 miRNAs in seven stages exceed 10,000 RPM, including chi-miR-1, chi-miR-206, chi-miR-148a-3p, chi-miR-381, chi-miR-127-3p, chi-let-7i-5p, chi-miR-26a-5p, chi-miR-10b-5p, chi-miR-378-3p, chi-let-7f-5p, chi-miR-99a-5p, chi-miR-199a-3p, chi-let-7g-5p, chi-miR-133a-3p, chi-miR-379-5p, chi-miR-143-3p, chi-miR-532-5p, chi-miR-411a-5p, and chi-miR-126-3p (Supplementary Table S3 and Figures 4A,B). And miR-1 had the highest expression in all stages, accounting for 23.36% of all expressions (Supplementary Table S3). Their constant expression suggested that they may perform essential cellular biological functions during skeletal muscle development to maintain cell integrity. Besides, we analyzed DEmiRNAs that were uniquely expressed in successive developmental stages to determine the unique role of DEmiRNAs at the seven stages ( Figure 4C). Specifically, 80 DEmiRNAs were found only in B1 and B90, precluding the prenatal DEmiRNAs (Supplementary Table S4 and Figure 4C). Among these DEmiRNAs, 34 and 46 were upregulated and downregulated, respectively. The upregulated DEmiRNAs were mainly involved in "HIF-1 signaling pathway, " "ErbB signaling pathway, " and other pathways that are associated with skeletal muscle development. Conversely, downregulated DEmiRNAs were involved in the "regulation of actin cytoskeleton, " "fatty acid biosynthesis, " and others pathways, which were related to intramuscular fat development (Supplementary Table S5). The earliest comparison group (F45 vs F65) participated in the "glycerophospholipid metabolism, " "GnRH signaling pathway, " and "citrate cycle (TCA cycle)" (Supplementary Table S5). Likewise, 23 DEmiRNAs were only found in the F90 and F120 comparison groups and were involved in "ribosome, " "Ras The matrix in the middle represents different sets. There is only one black dot in the same column, indicating the set expressed only in the group (left). The same column has two black dots, indicating the intersection of these two groups. Three black dots in the same column represent the intersection of these three groups and so on. The histogram on the right is the number of DEmiRNAs in the corresponding set. signaling pathway, " "rap1 signaling pathway, " and "glutathione metabolism" (Supplementary Tables S4, S5 and Figure 4C). And 28 DEmiRNAs only expressed between F120 and F135, which were involved in "SNARE interactions in vesicular transport, " "hedgehog signaling pathway, " "glycosaminoglycan biosynthesis-heparan sulfate/heparin, " and "lysosome" (Supplementary Tables S4, S5 and Figure 4C). In addition, eight DEmiRNAs in the F135-vs-B1 comparison group participated in the most signaling pathways, among which are "glycine, serine and threonine metabolism, " "arginine and proline metabolism" associated with amino acid metabolism, and "fatty acid degradation, " "fat digestion and absorption, " and "glycerolipid metabolism" associated with intermuscular fat production (Supplementary Tables S4, S5 and Figure 4C). These DEmiRNAs may represent the initiation/termination of physiological processes at a particular developmental stage in the skeletal muscle.

Trend Analysis Reflects the Role of miRNAs in Goat Skeletal Muscle
In order to understand the dynamic changes of DEmiRNA levels during goat skeletal muscle development more clearly, we selected DEmiRNAs with P-adj < 0.05 and |fold change| > 2 for trend analysis. A total of 305 DEmiRNAs were clustered into 12 expression patterns, three of which were significantly enriched, including profiles 16, 11, and 19 (Supplementary Figure S1). Profile 16 contained the most DEmiRNAs (188) and had a high expression level before birth (Figures 5A,B). They were enriched (P < 0.05) in 242 GO terms that mainly related to various enzymes, such as "SUMO-specific protease activity, " "SUMO binding, " and "isopeptidase activity" (Figure 5C and Supplementary Table S6). And they were involved in six enriched (P < 0.05) pathways, including "ribosome, " "mTOR signaling pathway, " "endocrine and other factor-regulated calcium reabsorption, " and other pathways ( Figure 5D and Supplementary Table S7). Among the DEmiRNAs in profile 16, 126 DEmiRNAs targeted 2,903 mRNAs. Five DEmiRNAs with the most target genes were identified as key DEmiRNAs for this profile, namely, novel-283, chi-miR-107-3p, chi-miR-328-3p, chi-miR-324-3p, and chi-miR-432-3p ( Figure 5E). Profile 11 (81 DEmiRNAs), which gradually decreased in all stages, targeted 1,340 mRNAs (Figures 6A,B). Functional analysis revealed that these target genes were enriched in 110 terms, 54 of which were related to response, enzyme activity, and the metabolic process, including "cellulase activity, " "cellular response to hormone stimulus, " and "cellular response to steroid hormone stimulus" (Figure 6C and Supplementary Table S8). Besides, some of the enriched functional terms of the skeletal muscle associated with proliferation, such as "cell proliferation, " "intermediate filament, " and "intermediate filament cytoskeleton" (Supplementary Table S8). And pathway analysis was enriched in "ribosome, " "insulin secretion, " "phototransduction, " and "sulfur relay system" (Figure 6D and Supplementary Table S9). The miRNA-mRNA network of profile 11 showed five key miRNAs, namely, chi-miR-665, chi-miR-412-3p, chi-miR-1271-3p, chi-miR-1306-3p, and chi-miR-767 ( Figure 6E). Furthermore, profile 19 contained 36 DEmiRNAs with 301 targeting mRNAs, and their expression levels rose in all stages (Figures 7A,B). The 58 functional enriched terms are related to material transport and structural stability, including "amino acid transmembrane transport, " "maintenance of protein location, " and "cytoskeleton The top 10 KEGG enrichment terms of profile 16. The bubble size represents the number of genes contained. P-value is from 0 to 0.06; red represents low P-value, and green represents high P-value. The yellow box represents the pathway with a P-value < 0.05. (E) miRNA-mRNA network for profile 16. The red triangles represent the DEmiRNAs, the yellow circles represent mRNAs, and the size represents the number of miRNA target genes. The larger the size, the greater the number of miRNA target genes. The blue word represents the top five nodes and is also the key miRNAs.

QRT-PCR Validation of DEmiRNAs
To validate data from RNA-seq, we validated five DEmiRNAs using qRT-PCR, including two miRNAs with high expression and one key miRNA in each expression profile. All of the verified miRNA expression levels showed a trend consistent with the amount of RNA-seq expression (Figure 8). miR-1 maintained a high expression level after F135, while miR-381 had a higher expression level from F120 to B1 (Figures 8A,B). In addition, chi-miR-328-3p, chi-miR-767, and chi-miR-150 were the key miRNAs of profiles 16, 11, and 19, respectively, and their qRT-PCR expression and RNA-seq expression were also consistent (Figures 8C-E).

DISCUSSION
As a branch of non-coding RNAs, miRNAs actively participate in the regulation of skeletal muscle development. With the development of sequencing technology, it is found that miRNAs in animal transcriptomes are more abundant than previously thought. There are now thousands of such molecules in skeletal muscles of many species, such as humans (Horak et al., 2016), sheep (Liu et al., 2019), and goat (Guo et al., 2016). Horak et al. (2016) reviewed several representative miRNAs found in basic research and sequencing that were specifically expressed in skeletal muscle and played an important role in proliferation, differentiation, and regeneration, including miR-1, miR-133a, miR-133b, miR-206, miR-208b, miR-486, and miR-499. These miRNAs are DEmiRNAs in this study. Liu et al. (2019) analyzed the GO and KEGG enrichment terms for the 345 DEmiRNAs of sheep's longissimus dorsi between fetuses and adults and showed the regulatory network of 29 DEmiRNAs and 298 target genes. In addition, Wang et al. (2014) revealed the expression profiles of goat miRNAs in the fetal and 6-month-old stages, identified 336 DEmiRNAs, and found that miR-424-5p and miR-29a may regulate muscle development. Furthermore, Guo et al. (2016) explored the miRNA expression profiles and their characteristics in four stages of goat skeletal muscle and the functions and pathways of DEmiRNA (RPM > 100) target genes at each stage of transformation. In addition, the relative weight gain rate of goats was accelerated before 90 days, while the absolute weight gain was accelerated after 90 days, especially after 120 days (Zhao et al., 2012). And miRNAs may change in the skeletal muscle before and after birth. Therefore, we chose seven-stage AWGs, namely, F45, F65, F90, The top 10 KEGG enrichment terms of profile 11. The bubble size represents the number of genes contained. P-value is from 0 to 0.06; red represents low P-value, and green represents high P-value. The yellow box represents the pathway with a P-value < 0.05. (E) miRNA-mRNA network for profile 11. The interpretation of the miRNA-mRNA network is the same as in Figure 5. F120, F135, B1, and B90. In our study, we not only analyzed the expression profile of miRNAs during the seven skeletal muscle developmental stages before and after goat birth and explored the possible role of miRNA target genes in each transformation stages and clustered trend expression patterns of DEmiRNAs but also performed functional and pathway analyses of enrichment trends to better understand the possible role of miRNAs with the same trend throughout development and screened out key miRNAs in each trend based on the miRNA-mRNA network.
All identified miRNAs (421 known miRNAs and 228 goat novel miRNAs) were concentrated between 21 and 23 nt, and bases were at most U, which was consistent with reported characteristics of miRNAs (Bartel, 2004). From this, we verified the reliability of the data and then analyzed the differences. A total of 420 DEmiRNAs were found in any two comparison groups. Among them, 19 DEmiRNAs expressed more than 10,000 RPM. Many of these miRNAs were involved in muscle development. miRNA-1 and miRNA-133a served as examples; they mediated the Dlk1-Dio3 Mega gene cluster inhibition to achieve metabolic maturation during muscle stem cell differentiation (Wust et al., 2018). miRNA-1 and miRNA-133a-3p were expressed at a low level in the early stage and an extremely high level after F120, indicating that they were essential factors regulating muscle differentiation and maturation after F120. And miR-206 promoted myoblast differentiation and muscle fiber hypertrophy, and it played an active role in the treatment of skeletal muscle injury (Ma et al., 2015). miR-206 was also highly expressed after F120, indicating that the progresses of differentiation and hypertrophy of goat skeletal muscle increased at F120. In addition, miR-148a-3p regulated proliferation and apoptosis of bovine muscle cells by targeting KLF6 (Song et al., 2019). It had a high expression level in the early stage, indicating that it might promote muscle cell proliferation and inhibited the process of apoptosis before goat F120. These results suggested that miRNAs that highly expressed in F45 to B90 may have an essential effect on skeletal muscle development.
To determine the unique role of DEmiRNAs in the seven developmental stages of skeletal muscles, we analyzed the unique DEmiRNAs in two consecutive stages of comparison groups. The most diverse DEmiRNAs in the B1 and B90 comparison groups are upregulated in the "HIF-1 signaling pathway, " "ErbB signaling pathway, " and others. The ErbB signaling pathway was thought to play an important role in myocyte proliferation, and the HIF-1 signaling pathway improved cell hypertrophy The bubble size represents the number of genes contained. P-value is from 0 to 0.05; red represents low P-value, and green represents high P-value. The yellow box represents the pathway with a P-value < 0.05. (E) miRNA-mRNA network for profile 19. The interpretation of the miRNA-mRNA network is the same as in Figure 5.  (Lebrasseur et al., 2003;Tan et al., 2010;Shang et al., 2013). This suggested that upregulated miRNAs in the B1 and B90 comparison groups that were uniquely expressed after birth released hypertrophic inhibition and regulated myocyte proliferation. The "TCA cycle" of the earliest comparison group (F45 vs F65) helped maintain muscle energy requirements and tissue homeostasis (Capitanio et al., 2017). The Ras signaling pathway involved in F90 and F120 comparison groups also played an important role in skeletal muscle development. The dysregulated Ras/MAPK pathway activity directly inhibited skeletal muscle development (Tidyman et al., 2011). The F135-vs-B1 comparison group was involved in a variety of amino acid signaling pathways, including glycine, serine, and threonine metabolism, all of which played an important role in skeletal muscle cells (Ost et al., 2015;He et al., 2017;Chen et al., 2018). For example, serine/threonine kinase 40 (Stk40) acted as a positive regulator of skeletal myoblast differentiation in culture and fetal skeletal muscle formation (He et al., 2017).
Finally, we clustered all miRNAs and found three major trends, namely, profiles 11, 16, and 19. For profile 16, they have a high level of expression in the prenatal phase. The enrichment of GO and KEGG may mean that prenatal miRNAs have reduced inhibition of skeletal muscle hypertrophy and differentiation-related pathways, and the mRNA expression associated with these pathways begins to rise. chi-miR-107-3p had a high expression level (RPM > 100) before F120, which can target 173 mRNAs. For example, its target gene DLL-1 interacts with the NOTCH receptor to reduce the transmission of Notch signaling, leading to skeletal muscle hypertrophy (Al Jaam et al., 2016). Another target gene, DOCK3, that overexpressed it could regulate PTEN/AKT signaling, thereby regulating muscle hypertrophy and growth and inducing apoptosis (Alexander et al., 2014). chi-miR-324-3p, another key miRNA of profile 16, targeted AKT1, which promoted type IIb muscle growth (Cheng et al., 2015). And AKT1 also promoted load-induced muscle hypertrophy by regulating satellite cell proliferation (Moriya and Miyazaki, 2018). Besides, inhibition of chi-miR-328-3p targeting gene PLD prevented mechanically induced PA elevation and activation of mTOR signaling (Hornberger et al., 2006). Moreover, novel-283 targeted AMOT and BCL9. Amot was a necessary factor to regulate the expansion and self-renewal of muscle stem cells (Li and Fan, 2017). BCL9, as an essential component of canonical Wnt signaling, mediated the differentiation of myogenic progenitors during muscle regeneration (Brack et al., 2009). And chi-miR-432-3p targeted NPRL2, and its reduced expression level would make muscle fibers larger (Dutchak et al., 2018). These late-decreasing miRNAs played essential roles in muscle hypertrophy. Key DEmiRNAs were also screened in profile 11, which was declining from F45. Among them, the targeted gene ANGPT4 of chi-miR-767, the targeted gene POLR3H of chi-miR-665, and the targeted gene TEPSIN of chi-miR-1271-3p were all involved in the enriched functional term "cell differentiation." Targeted gene MARK2 of chi-miR-1306-3p regulated muscle stem cell polarity to ensure asymmetric separation (Dumont et al., 2015). Furthermore, five key DEmiRNAs were also selected in the rising spectrum at all stages. The targeted genes of these key DEmiRNAs were involved in the "cytoskeleton organization, " including the targeted genes GDF7, LOC102169433, and VSIG10 of chi-miR-150; ANKRD13D and MLPH of chi-miR-193b-5p; KCNMA1 and LRRC73 of chi-miR-361-3p; and C11H9orf16 of chi-miR-193a and chi-miR-193b-3p. They were also closely related to skeletal muscle development. For example, FGFR1, which was targeted by chi-miR-150, was important for FGF-mediated proliferation of skeletal muscle satellite cells and contributed to mitogenic effects (Yablonka-Reuveni et al., 2015). Moreover, SLIT1, also targeted by chi-miR-150, drove directional migration and differentiation of Robo2-expressing pioneer myoblasts (Halperin-Barlev and Kalcheim, 2011). In addition, the DEmiRNAs of the three expression profiles were most abundant in the ribosome, whether rising or falling. These data indicated that miRNAs in skeletal muscle development have different effects on ribosomes at different stages.

CONCLUSION
In this study, we identified miRNAs in goat skeletal muscle at five fetal stages and two kid stages by RNA-seq and systematically explored temporal trend expression profiles of skeletal muscle development. A total of 421 known miRNAs and 228 goat novel miRNAs were identified in the data. Moreover, 420 DEmiRNAs (364 known miRNAs and 56 goat novel miRNAs) were found in all comparison groups in the seven stages, and the average abundance of 19 miRNAs exceeded 10,000 RPM. The sustained high expression of these miRNAs indicates that they play important cellular biological functions during skeletal muscle development. In addition, DEmiRNAs were clustered into 12 trends and were mainly enriched in three temporal trend expression profiles (profiles 16, 11, and 19). Functional and pathway analyses indicated that the biological pathways regulated by miRNA expression profiles were directly related to temporal changes in goat skeletal muscle development. This study illustrated the role of miRNAs in different goat skeletal muscle stages, which provided valuable information for a better understanding of mammalian skeletal muscle development.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the RNAseq data was available in the SRA database under the accession number PRJNA553597.

ETHICS STATEMENT
This study was carried out following the principles of the Basel Declaration and recommendations of the Guide for the Care and Use of Laboratory Animals (http://grants1.nih. gov/grants/olaw/references/phspol.htm), the ethics committee of Anhui Agricultural University. Anhui Agricultural University Ethics Committee has approved the agreement under permit No. AHAU-AE2017-07.

AUTHOR CONTRIBUTIONS
YHL: conceptualization, funding acquisition, methodology, and supervision. YHL and QZ: data curation and formal analysis. QZ and LZ: investigation. MS and LZ: resources. QZ, MS, and LZ: validation. QZ and MS: visualization. QZ: writing -original draft. YHL, YZ, YaL, FF, YM, and XZ: writing -review and editing. All authors contributed to the article and approved the submitted version.  (201808340031), and Anhui University Talent Project (gxfx2017004). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2020.00684/full#supplementary-material FIGURE S1 | All trend clustering of DEmiRNAs. (A) We used the STEM to classify these genes into 12 categories. Color denotes enrichment.
TABLE S1 | The summary information of small RNA.