The experimental procedures were approved by the Institutional Animal Care and Use Committee at Shihezi University.

We collected the longissimus dorsi of three 1-year-old healthy female Kazak sheep were from a slaughterhouse. Surgery was carried out in a sterile environment to obtain the longissimus dorsi of three prenatal sheep. Tissue samples were immediately frozen in liquid nitrogen and stored at −80°C, respectively. The total RNA was subsequently extracted from the pooled samples using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions.

The cDNA was synthesized by reverse transcriptase polymerase chain reaction (PCR) and the cDNA library was sequenced using Hi Seq 2000 (Illumina, San Diego, CA, USA). To improve the quality of the sequencing, it was necessary to strictly filter reads that were contaminated with primer or adapter sequence or was of low quality. The small RNAs were mapped to the reference genome sequences of sheep by Bowtie software and their expression and distribution on the reference sequences were analyzed. Integrated with miRBase database, the sRNAs aligned to the reference sequences were found and separated as much as possible from other small RNA fragments.

Hairpin structure of the mRNA precursor can be used to predict new miRNAs. We integrated miREvo [7] and mirdeep2 software [8] to predict novel miRNA, according to the signature hairpin structure of the miRNA precursor.

To compare the differential expression of miRNAs in two libraries, we firstly normalized expression levels of miRNAs by using transcript per million (TPM) [9]. For biologically repeated samples, the analysis was performed using DESeq2 based on the negative binomial distribution [10]. For non-biologically repeated samples, TPM was used to standardize the read count data, followed by differential expression analysis with DEGseq [11]. The screening criteria for differentially expressed miRNAs were usually set by default to: p-value <0.05.

In this paper, the miRanda [12] was performed to predict targeted genes of differentially expressed miRNAs in skeletal muscle of prenatal and postnatal sheep. GOseq based on Wallenius non-central hyper-geometric distribution that eliminates gene length preference was used for gene ontology (GO) enrichment analysis [13]. And we performed Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of candidate targeted genes using KOBAS [14] software.

Five miRNAs (miR-181a, mir-199, mir-127, mir-133, and mir-1) were selected from differentially expressed miRNAs that were negatively correlated with their targeted genes and their targeted genes (calsequestrin 1 [CASQ1], muscleblind-like splicing regulator 1 [MBNL1], carbonic anhydrase 4 [CA4], insulin-like growth factor 2 [IGF2], and myelin expression factor 2 [MYEF2]) were characterized and quantified by quantitative PCR (qPCR). The cDNA was synthesized by reverse transcription for further experiments. We performed qPCR with each sample repeated three times. Quantitative real-time (RT)-PCR was performed as follows: 1 ng of cDNA was added to reaction mixture containing 5 μL 2×UltraSYBR Mixture (CWBIO, Beijing, China), 0.6 μL primers and water added to final volume of 10 μL. The miRNA detection program was 95°C for 10 min, followed by 45 cycles of 95°C for 25 s, 58°C for 10 s, 72°C for 10 s. The β-actin gene was used as reference gene. Genes expression levels were evaluated as 2^−ΔΔct. The primers used in this experiment were shown in Supplementary Table S1.

To detect miRNAs expression profile in sheep skeletal muscle, we constructed two sRNA libraries of the longissimus dorsi of sheep by high-throughput sequencing. We obtained about 16 and 15 million raw sequence reads per sample from prenatal library (Longissimus dorsi embryo [LDE, prenatal stage]) and postnatal library (Longissimus dorsi adult [LDA, pastnatal stage]), respectively (Table 1). After removing an adaptor and low quality reads, each library finally obtained 15.9 and 14.9 million clean reads with lengths between 18 to 35nt and length distribution peak at 22nt (Supplementary Figure S1), of which 88.81% and 94.24% were mapped to the sheep reference sequence, respectively. Among these sRNAs, most of them were aligned with non-miRNAs databases, such as rRNA, tRNA, snRNA (Supplementary Figure S2). Integrated with miRBase database, 1,086 of miRNAs were aligned to mature known miRNAs, which were first base biased U+G (Figure 1A). The remaining candidate sequences were used to predict potentially new miRNAs.

To further validate the miRNA sequences, new miRNAs were predicted using miREvo and mirdeep2 software. In result, a total of 40 novel RNAs were found in LDE and LDA (Supplementary Table S2). It was worth noting that the expression of Novel_1 was the highest in two different periods (Figure 2). Interestingly, through statistical analysis of the first base preference at different lengths of miRNAs, the data indicated that the first base was predominantly preferred for U, followed by G and C at two developmental stages (Figure 1B).

Statistical analysis of the expression levels of known and new miRNAs in each sample was performed and normalized with TPM. Among these 1,000 miRNAs, 615 were co-expressed, while 300 and 18 were expressed in LDE and LDA, respectively (Figure 3A). Correlations of miRNAs between LDE and LDA were lower (R² = 0.627, Figure 3B). To generate miRNAs with differential expressions during sheep muscle development, we detected the level of miRNAs expressions at two development stages and identified 254 up-regulated, 91 down-regulated (Figure 3C), and the top of 50 miRNA, significantly differentially expressed are shown in Figure 4. Meanwhile, hierarchical clustering of differentially expressed genes was performed, resulting in 16 subgroups. Among those subgroups, the largest subclass contained 50 differentially expressed genes (Figure 3D).

Targeted genes prediction of miRNAs was performed by the miRanda software, resulting in a predicted total of 2,833 targeted genes (Supplementary Table S1). We found that different miRNAs target the same gene, conversely, the same miRNA target different genes (Table 2). The distribution of differential expressed genes in GO can be used to study the physiological functions and the biological processes in different stages of sheep (Figure 5). In the skeletal muscle of sheep, the analysis for differentially expressed genes through GO enrichment analysis showed that 184 terms were preferentially enriched in cell components, including: endoplasmic reticulum, nuclear lumen, Golgi apparatus etc. while 204 terms were significantly enriched in molecular functions including: protein binding, nucleoside phosphate binding, adenyl ribonucleotide binding etc. And 982 terms were mostly enriched in biological processes, such as metabolic process, negative regulation of signaling, positive regulation of cellular catabolic process, enzyme linked receptor protein signaling pathway etc. These results suggested that miRNAs may be involved in the regulation of muscle development by participating in protein metabolic signaling pathways.

KEGG pathway analysis of miRNAs targets showed that 273 terms were enriched. Among these pathways, the main signaling transduction pathways associated with muscle development were: Wnt signaling pathway, transforming growth factor beta signaling pathway, Hippo signaling pathway, Thyroid hormone signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway etc. (Supplementary Table S3). These results indicated that miRNAs may be involved in regulation of the development and growth of skeletal muscle.

miRNAs, as regulatory factors, negatively regulate targeted gene expression at the post-transcriptional level by translational repression or by silencing its targeted genes. Therefore, we obtained some significant regulatory networks by analyzing down-regulation of miRNA versus up-regulated genes, and up-regulation of miRNA versus down-regulated genes, respectively. Through regulatory network analysis for down-regulated miRNA versus up-regulated mRNA, the results showed that 29 down-regulated miRNAs were correlated with 290 targeted genes (Figure 6A). Among them, 81% targeted genes may be regulated by miR-143, -miR-23a, miR-29b, miR-30c, miR-23b, miR-26a, miR-221, miR-150, miR-30a-3p, miR-191, miR-22-3p, miR-21, and miR-133. Among these targeted genes, oar-miR-133 might be regulating gene expression of insulin-like growth factor 2 (IGF2), intraflagellar transport 122, thyroid hormone receptor interactor 13, fat mass and obesity-associated protein (FTO), and ETS proto-oncogene 1 etc. In addition, novel_175 was correlated with insulin like growth factor 2 mRNA binding protein 1, fibroblast growth factor 14 (FGF14), v-myb avian myeloblastosis viral oncogene homolog-like 2, formin 2, meis homeobox 3 etc. Similarly, the network of up-regulated miRNA versus down-regulated mRNA in skeletal muscle is shown in Figure 6B. The result showed that 95 up-regulated miRNAs were associated with 62 targeted genes. MiR-106b may be regulate genes expression of paraoxonase 3, thioredoxin interacting protein, pyruvate dehydrogenase kinase 4, interferon related developmental regulator 1, muscular LMNA interacting protein, and tubulin alpha-8 chain. Novel_140 regulated carnitine O-acetyltransferase, solute carrier family 27 member 6, glycogenin 1, LIM and cysteine rich domains 1, solute carrier family 25 member 33, chloride voltage-gated channel 1, NADH:ubiquinone oxidoreductase core subunit S7, and myosin light chain 3.

To verify the reliability of the sequencing data, we performed qRT-PCR to detect the relative expression levels of five miRNAs selected from differentially expressed miRNAs that were negatively correlated with their targeted genes. The results showed that 3 miRNAs (miR-199a, miR-181a, and miR-127) were up-regulated, 2 miRNAs (miR-1 and miR-133a) were down-regulated, which were correlated with the RNA-seq analysis (Figure 7A). Moreover, we detected the expression of mRNA of these 5 miRNAs targeted genes by qPCR. There were 3 targeted genes (MBNL1, CASQ1, and CA4) that were down-regulated, and 2 targeted genes (MYEF2 and IGF2) that were up-regulated, which were negatively correlated with their corresponding miRNA (Figure 7B). Further, it has been proved from the results that the Illumina HiSeqTM2500 sequencing data was reliable.