Characterization of miRNA expression profiles
We calculated the length distribution of the total number of valid data, and certain RNA sequences, including rRNA, tRNA, snRNA, and snoRNA, were searched against and filtered from the raw reads by aligning the sequences to RFam and Repbase databases. The length distribution showed that over half of the clean reads (52.6%) were 22 nt in length, which was a good cue that most of the clean reads obtained were miRNA sequences. All clean reads were aligned to the sheep genome and mammalian miRbase using the ACGT-miR101 program. Finally, 4,752 miRNAs were detected in total, and 2,275 of them have never been reported.
Identification of DE-miRNAs and functional analysis
The DE-genes were screened with p ≤ 0.05 as a threshold. In the expression profile, a total of 505 miRNAs were differentially expressed. According to the expression characteristics of profiles, there were five categories: incremental type (107), decreasing type (151), high-low-high type (91), and low-high-low type (45), and irregularity type(111). We performed a pairwise comparison of three groups on the basis of filter criteria(D85N vs D105N, D105N vs D135N, and D85N vs D135N). There were 63 up-regulated miRNAs and 43 down-regulated miRNAs in D85N vs D105N, 106 up-regulated miRNAs and 123 down-regulated miRNAs from D105N to D135N, and 172 up-regulated miRNAs and 112 down-regulated miRNAs between D85N and D135N. We then performed overlapping statistics on these data, and 16 of the DE-miRNAs were the same in three groups (Fig. 1). The detailed information of these 16 miRNAs was extracted from the miRNA profile (Table 1). According to the sequence, it is found that 2 of them are the same miRNA. Duplicate data was deleted and 14 miRNAs were used for subsequent network construction.
Table 1
The detailed sequencing information of 16 miRNAs
miR_name | miR_seq | Expression level | Pvalue | Regulation Type |
hsa-miR-410-5p | AGGTTGTCTGTGATGAGTTCG | middle | 1.01E-08 | up |
bta-miR-378c_R-1_1ss20GT | ACTGGACTTGGAGTCAGAAT | middle | 1.62E-06 | down |
oar-miR-150 | TCTCCCAACCCTTGTACCAGTG | middle | 4.40E-06 | down |
bta-miR-193a-5p_R-3 | TGGGTCTTTGCGGGCGAGA | middle | 4.92E-06 | down |
oar-miR-3959-3p | TGTATGTCAACTGATCCACAGT | middle | 6.13E-06 | up |
hsa-miR-378d | ACTGGACTTGGAGTCAGAAA | middle | 4.73E-05 | down |
oar-miR-23b_R + 6 | ATCACATTGCCAGGGATTACCACT | high | 5.53E-05 | down |
chi-miR-133a-5p | AGCTGGTAAAATGGAACCAAAT | high | 1.68E-04 | down |
chi-miR-365-3p | TAATGCCCCTAAAAATCCTTAT | high | 1.70E-04 | down |
bta-miR-365-3p | TAATGCCCCTAAAAATCCTTAT | high | 1.70E-04 | down |
oar-miR-493-3p | TGAAGGTCTACTGTGTGCCAGG | high | 1.71E-04 | up |
hsa-miR-135a-3p_L + 1R-1 | ATATAGGGATTGGAGCCGTGGC | middle | 2.84E-04 | up |
oar-miR-134-5p_R + 1 | TGTGACTGGTTGACCAGAGGGT | middle | 2.92E-04 | up |
chi-miR-450-5p_R + 2 | TTTTGCGATGTGTTCCTAATAT | high | 5.63E-04 | up |
bta-miR-450a | TTTTGCGATGTGTTCCTAATAT | high | 5.63E-04 | up |
bta-mir-2284z-2-p3_1ss22GT | AAAAACCCAGATGAACTTTTTT | middle | 1.16E-03 | up |
To better understand the functions of DE-miRNAs, we conducted GO term and KEGG pathway analyses. A total of 643 unique biological processes and 73 enriched pathways were identified (p ≤ 0.05), including signal transduction processes such as MAPK, Wnt, regulation of actin cytoskeleton, Ras, focal adhesion and the ErbB signaling pathway. Based on the top 20 GO enrichment analysis of miRNAs (Fig. 2a), including biological process, cellular component, molecular function. In biological process, they were riched in protein ubiquitination involved in ubiquitin-dependent protein catabolic process, protein dephosphorylation, protein glycosylation, microtubule cytoskeleton organization, regulation of gastrulation, peptidyl-tyrosine dephosphorylation and actin filament organization, which suggested that spatio-temporal embryonic development was closely related to intracellular protein breakdown and maintenance, energy and material metabolism, signal transduction, and et al. We then further analyzed the 16 miRNAs selected above to better understand their functions. In the top 20 KEGG pathway(Fig. 2b), the most highly enriched pathway was nucleotide excision repair. Furthermore, many pathways related to metabolism were also enriched, such as Wnt signaling pathway, insulin signaling pathway, focal adhesion, and ErbB signaling pathway.
miRNA-mRNA networks
With 505 differential expressed miRNAs, more than 100,000 target genes were screened in the database. Therefore, we focused on the 14 important miRNAs selected above for target prediction and screened 1,137 target genes in total. Then, integrating the differential expressed mRNAs in profile, 379 target genes were retained, and 944 miRNA-mRNA interaction pairs related to skeletal muscle fiber development were constructed (Fig. 3). Each miRNA regulates multiple mRNAs and vice versa. A total of 253 target genes were obtained from 7 down-regulated miRNAs and a total of 126 target genes were obtained from 7 up-regulated miRNAs. Among them, oar-miR-150 regulates the most target genes, with a total of 100. The average number of target genes regulated by up-regulated and down-regulated miRNAs is more than 35, except for bta-miR-450a.
miRNA-TG-pathway network
On the basis of miRNA-mRNA analysis, we performed an analysis of the miRNA-TG-pathway network in order to find the biological processes and signaling pathways that the miRNA target genes are mainly involved in. From the results of BINGO, it can be seen that the DE-genes are mainly concentrated in 198 biological processes (Additional file 2), which mainly cover 5 aspects, such as multicellular organismal process, intracellular signal transduction, reproductive process, animal organ development and metabolic process, of which muscle tissue development, actin filament-based process, muscle structure development, and muscle cell differentiation processes have attracted our attention (Fig. 4a). In the diagram drawn by ClueGo and CluePedia (Fig. 4b), 5 pathways were enriched, cytoskeleton organization, positive regulation of synaptic transmission, regulation of small GTPase mediated signal transduction, peptidyl-serine phosphorylation, regulation of neuron differentiation, among which More than 47% of genes are involved in the cytoskeleton tissue process. Finally, the DE-genes involved in these five pathways were extracted, and miRNAs from these DE-genes were extracted to construct a miRNA-TG-pathway regulatory network. In the miRNA-TG-pathway network (Fig. 5), the enriched pathways are Insulin signaling pathway, Wnt signaling pathway, Cell cycle, Focal adhesion, Adherens junction, Ubiquitin mediated proteolysis, ECM-receptor interaction, ErbB signaling pathway, Tight junction, p53 signaling pathway. Among them, Focal adhesion pathway enriched the most differentially expressed target genes. Seven miRNAs enriched in this network were used for in-depth analysis.
Integral lncRNA-miRNA-mRNA interaction networks
Integrating the current miRNA data and the previously obtained lncRNA data,we performed ceRNA analysis. The differentially expressed lncRNAs were screened to predict their target miRNAs, and the data containing the seven miRNAs selected from miRNA-TG-pathway network were extracted for the construction of lncRNA-miRNA-mRNA networks. Among them, 14 differentially expressed lncRNAs were able to adsorb oar-miR-493-3p, bta-miR-450a, oar-miR-23b_R + 6, hsa-miR-410-5p and oar-miR-3959-3p. However, chi-miR-365-3p and chi-miR-133a-5p were not able to bind to any differentially expressed lncRNAs. Then, 1,598 edges in lncRNA-miRNA-mRNA networks were constructed based on the ceRNA mechanism(Fig. 6a). By analyzing its topological characteristics, the top 25 nodes with the highest degrees were identified (Fig. 6b). Among them, oar-miR-493-3p has the highest frequency in the network, and is related to genes (TEAD1, MAP1B, MCM9, BRI3BP, MPP7, DGKE, MAN1A2 and ZBTB34) and lncRNA (MSTRG.4058, MSTRG.4324, MSTRG.2252, MSTRG. 3879 and MSTRG.1470), which is considered to be of great significance in ceRNA networks. Furthermore, oar-miR-3959-3p and hsa-miR-410-5p also regulate multiple genes and lncRNA simultaneously, so oar-miR-493-3p, oar-miR-3959-3p and hsa-miR-410-5p were selected as important miRNAs. Similarly, according to the degree value of lncRNAs in the network, MSTRG.3533, MSTRG.4324 and MSTRG.1470 are selected as important candidate lncRNAs.
miRNA-TG-TF network in WGCNA modules
In this section, we combine the previous WGCNA results[14] to perform in-depth analysis of the screened miRNAs and target genes above. The WGCNA analysis of lncRNA and mRNA was clustered into 25 modules.It can be seen from the that of all the modules, the turquoise module has the highest frequency in the previously study, so it has the most enrichment significance. Most differentially expressed RNAs in the ceRNA network fall into this turquoise module. The functional enrichment analysis of DE-genes in the turquoise module shows that these genes participate in many signaling pathways, including ECM-receptor interaction, Adhesion, mTOR signaling pathway, TypeⅡdiabetes. It is consistent with the above results. The hub genes were extracted from it and predict their TF using USCS. Then, integrating the interaction pairs, a miRNA-TF-gene network was constructed(Fig. 7a). The top 30 nodes with the highest degrees were identified based on its topological characteristics and were constructed a sub-network (Fig. 7b). It can be seen that 4 genes(TEAD1, ZBTB34, GSK3B and POGLUT1) and 2 miRNAs(oar-miR-493-3p and oar-miR-3959-3p) play important roles in sub-network. GSK3B and POGLUT1 as DE-genes that act as TFs are shown to be targeted by oar-miR-3959-3p and seen to regulate the expression of many TFs (i.e., C / EBPbeta, Sp1, NF1, TFIID and PR B). TEAD1 and ZBTB34 act as TFs are shown to be targeted by oar-miR-493-3p, they can also regulate TFs. In addition, C / EBPbeta, Sp1, NF1, TFIID and PR B can interact with at least three DE-genes. Therefore, according to their functions in skeletal muscle development, C / EBPbeta, TFIID and PR B were selected as important muscle development-related transcription factors.
Validation of RNA-Seq data results
The validation results for the five miRNAs selected to substantiate the accuracy of sequencing are displayed in Fig. 8. The results indicate that there is a similar expression pattern of miRNAs generated from RNA-Seq and qRT-PCR data.