miRNA sequencing and quality assessment results
All sequencing data for each small RNA library are presented in Additional file 1: Table S1. We mapped these clean reads using the sheep genome from Ensembl to further elucidate the possible mechanisms underlying the diversity of small RNAs involved in HF development. Mapping the total reads to parts of the genome resulted in more than 80% of the clean reads in the 18 libraries (Additional file 1: Table S1), and the distribution of small RNA length differed among the eighteen libraries. We filtered raw data, and the results showed that the coverage of clean reads Q30 of miRNA was more than 75% (Additional file 1: Table S1). In all libraries, most of the clean reads were 22 nucleotides long (Additional file 2: Fig. 1b). The sequencing depth and coverage of samples met the needs of this experiment and could be used for subsequent experiments.
qRT–PCR Validation of miRNA-Seq
We compared the qRT–PCR results with the miRNA-seq results (Fig. 1). We found that the qRT–PCR expression levels of miR-148a in three periods (E135, D7 and D30) and miR-148a and of miR-154a-5p in four periods (E65 to E135) in skin were inconsistent with the sequencing results. The qRT–PCR expression levels of the remaining seven miRNAs were consistent with the miRNA-seq results. Therefore, we concluded that the sequencing results were accurate and reliable.
Identification of differentially expressed miRNAs
A total of 87 DE-miRNAs and 446 novel DE-miRNAs were identified in the six HF development stages of SBM sheep. Previous research by the group identified 7879 DE-mRNAs [53]. A total of 11 DE-miRNAs and 104 novel DE-miRNAs were identified in E65 vs. E85, of which 5 DE-miRNAs were upregulated and 6 were downregulated, while 39 novel DE-miRNAs were upregulated and 65 were downregulated (Fig. 2a). When E85 was compared to E105, 15 DE-miRNAs (including 11 upregulated and 4 downregulated DE-miRNAs) and 104 novel DE-miRNAs (including 62 upregulated and 42 downregulated novel DE-miRNAs) were detected. In total, 43 DE-miRNAs (6 upregulated and 37 downregulated) and 97 novel DE-miRNAs (50 upregulated and 47 downregulated) were found between E105 and E135. Between E135 and D7, 9 DE-miRNAs (5 upregulated and 4 downregulated) and 105 novel DE-miRNAs (78 upregulated and 27 downregulated) were identified. When D7 was compared to D30, 21 DE-miRNAs (1 upregulated and 20 downregulated) and 48 novel DE-miRNAs (21 upregulated and 27 downregulated) were identified. E65 vs. D30 had the most DE-miRNAs, with 60, and E135 vs. D7 had the fewest DE-miRNAs, with only 9. This result indicated considerable differences between the beginning of HF development and postnatal development. Venn diagram analysis (Fig. 2b) indicated that miR-29a was differentially expressed in the five groups (except E65 vs. E85). miR-23b, miR-655-3p, miR-431, miR-410-3p and miR-29b were differentially expressed in four different periods. Likewise, 44 DE-miRNAs were specifically expressed in only one period, of which E65 vs. E85 had 2, miR-181a and miR-433-5p, respectively. There were 3 miRNAs (miR-154b-3p, miR-323b and miR-655-5p) specifically expressed in E85 vs. E105. Eleven E105 vs. E135 (miR-136, miR-380-5p, miR-411b-5p, etc.) were specifically expressed. Two (miR-1197-5p and miR-410-5p) and one (487a-3p) miRNAs were specifically expressed in E135 vs. D7 and D7 vs. D30, respectively. E65vs. D30 is the most specifically expressed miRNA, with 25.
Enrichment analysis of miRNA target genes
To further elucidate the functions of the DE-miRNAs, we also examined the top 10 GO functional enrichment BP (biological processes) terms (Fig. 3, Additional file 3: Fig. S2) and KEGG pathways (Additional file 4: Fig. S3) of target genes in adjacent comparison groups. For E65 vs. E85, the GO terms were enriched mainly in the positive regulation of vacuole organization, innate immune response, extrinsic apoptotic signaling pathway and negative regulation of Notch signaling pathway. The GO terms E85 vs. E105 were enriched in cell proliferation involved in outflow tract morphogenesis and segmentation. For E105 vs. E135, the terms were mainly concentrated in regulation of apoptotic process and regulation of programmed cell death. For E135 vs. D7, the terms were mainly enriched in BPs related to negative regulation of the Notch signaling pathway, positive regulation of vacuole organization, and outflow tract septum morphogenesis. For D7 vs. D30, the terms were mainly enriched in cell cycle checkpoint, tube closure and chordate embryonic development. In addition, in the comparative analysis between 30 days after birth (D30) and the initial stage of HF development (E65), the target genes were enriched in outflow tract septum morphogenesis and chordate embryonic development. In summary, the target genes in E65 vs. E85 were directly related to cellular immunity and apoptosis. The target genes in E105 vs. E135 were significantly correlated with apoptosis. The target genes in E135 vs. D7 were related to kidney development, and those in E65 vs. D30 were related to cardiac and neurological development.
We performed KEGG enrichment analysis and found that the target genes were mainly enriched in the proteasome, AMPK signaling pathway, protein processing in endoplasmic reticulum and other pathways. Similarly, it was also enriched in the neomycin, kanamycin and gentamicin biosynthesis pathways (Additional file 4: Fig. S3).
DE-mRNAs and DE-miRNAs cluster analysis
To screen mRNAs and miRNAs, according to heatmap analysis, the DE-mRNAs at the first two time points (E65 and E85) and last four time points (E105, E135, D7 and D30) were combined as Stage A and Stage B, respectively (Fig. 4a). Similarly, the DE-miRNAs at the first three time points (E65, E85 and E105) and the later three time points (E135, D7 and D30) were combined as Stage A and Stage B (Fig. 4c). Then, we generated a Venn diagram for Stage A and Stage B. The number of Stage A DE-mRNAs and DE-miRNAs was 1562 and 21, respectively, and the number of Stage B DE-mRNAs and DE-miRNAs was 2100 and 28, respectively. The numbers of DE-mRNAs and DE-miRNAs in the two groups were 428 and 6. The DE-mRNAs that appeared in both groups were eliminated, and the remaining 1134 DE-mRNAs were identified as candidate genes associated with their respective HF development stage (Stage A). The remaining 1672 DE-mRNAs are candidate genes associated with the hair follicle maturation stage (Stage B) (Fig. 4b). The remaining 15 DE-miRNAs were candidates associated with stage A. The remaining 22 DE-miRNAs were candidates associated with stage B (Fig. 4d).
Enrichment analysis of mRNA Stage A and Stage B
To further study the biological functions of DE-mRNAs related to the development of follicles in SBM during the fetal period, GO functional analysis and KEGG pathway analysis were performed on stage A DE-mRNAs (Fig. 5). In biological processes, the DE-mRNAs were significantly enriched in negative regulation of the canonical Wnt signaling pathway, HF development, establishment of the skin barrier, skin development, cellular response to transforming growth factor beta stimulus, negative regulation of the BMP signaling pathway, negative regulation of epithelial cell proliferation, epithelial-to-mesenchymal transition, positive regulation of epidermal cell differentiation and HF morphogenesis (Fig. 5a). These biological processes have been reported to be associated with HF morphogenesis, skin development or even hair formation. We visualized the genes involved in these biological processes and found many DE-mRNAs. However, we focused on the 30 most meaningful genes (e.g., TGFβ2, NOTCH1, WNT10A, ACVR1B, SOSOTCD1, WNT4, and BMP2) (Fig. 5b). It is presumed that the initiation of HF development is mainly driven by the above factors.
A list of significant pathways was obtained by KEGG analysis of the 1134 DE-mRNAs in Stage A (Fig. 5c). Four significant pathways were identified: the Hippo signaling pathway, PI3K-Akt signaling pathway, TGF-β signaling pathway and Hedgehog signaling pathway. These pathways have also been reported to be significant in HF development in previous studies. We further analyzed the network of KEGG pathways related to HF development in Stage A and found 20 genes, including TGFβ2, WNT10A, PTCH1, BMP2, ACVR1B, WNT5A, WNT16, LAMA5, WNT10A, and WNT4, that may be involved in HF development. In summary, we found that 19 genes (LAMA5, ITGB4, ITGA3, THBS1, NOG, ACVR1B, INHBA, TGFβ2, BMP2, BMP5, WNT10A, SOX2, WNT16, WNT14, WNT5A, PTCH1, PTCH2, SHH, and FZD1) play a key role in Stage A of HF development.
A list of significant pathways was obtained by KEGG analysis of 1672 DEGs in Stage B (Fig. 6). We were surprised to find that only one biological process reported to be associated with HF development was enriched in Stage B, positive regulation of NF-kappa B transcription factor activity (Fig. 6a). Among the identified pathways, the PPAR signaling pathway has been most studied in the development of HF (Fig. 6b). Our study suggested that the genes in this pathway are associated with the late stage of HF development and maturation.
Networks of DE miRNA-DE mRNA and miRNA target DE-mRNA analysis
Predicted target genes were combined with transcriptome profiling data to construct an mRNA–miRNA network associated with sheep HF growth and development. Significant mRNA–miRNA pairs in Stage A were selected to identify a total of 17 target genes of the 9 known sheep DE-miRNAs. The mRNA–miRNA network is shown in Fig. 7a.
To understand the roles of mRNAs and miRNAs in controlling skin morphogenesis and HF development, the expression of miR-133 and miR-23b and their target genes ACVR1B and WNT10A was evaluated in the skin of sheep at different stages of embryonic development and at postnatal D7 and D30. A low level of miR-133 expression was detected in embryonic skin during the onset of HF development from E65 to E85; its expression was significantly increased in E105 skin and then decreased significantly by E135. ACVR1B exhibited the opposite expression pattern (Fig. 7b). miR-23b expression gradually increased and was maximally expressed on D30, whereas WNT10A expression showed the opposite pattern from E65 to E105 (Fig. 7c). Taken together, these findings suggest that miR-133 targets ACVR1B, while miR-23b targets WNT10A. To detect the targeting relationship between these miRNAs and genes, we constructed luciferase reporter vectors psiCHECK2-WNT10A and psiCHECK2-ACVR1B to transfect 293T cells.
Validation of miR-133 and miR-23b functions in sheep dermal fibroblasts
To validate ACVR1B and WNT10A as direct targets of miR-133 and miR-23b, respectively, ACVR1B and WNT10A luciferase reporters were constructed with wild-type and mutated ACVR1B and WNT10A 3’UTRs (Fig. 8a). These reporters were cotransfected with miR-133 and miR-23b mimics or a negative control (mimic-NC) into 293T cells. Luciferase reporter activity was measured in dual luciferase assays. As shown in Fig. 8, the luciferase activity was significantly decreased when cotransfected with miR-133 and miR-23b mimics compared to either mutant or empty controls. Overexpression of miR-133 and miR-23b mimics decreased the luciferase activity of psiCHECK2-ACVR1B (Fig. 8b) and psiCHECK2-WNT10A (Fig. 8c).
We verified the network related to ACVR1B and WNT10A and showed that ACVR1B and WNT10A are target genes of miR-133 and miR-23b, respectively. To test whether ACVR1B and WNT10A can be regulated by miR-133 and miR-23b, the mRNA levels of ACVR1B and WNT10A in sheep dermal fibroblasts transfected with miR-133 and miR-23b mimics and inhibitor were determined by qRT–PCR. The miR-133 and miR-23b mimics dramatically decreased the mRNA levels of ACVR1B and WNT10A, while miR-133 and miR-23b inhibitors increased the ACVR1B (Fig. 8d) and WNT10A (Fig. 8e) mRNA levels in sheep dermal fibroblasts. These results were consistent with the observed decrease in WNT10A expression (Fig. 8c). To further test whether ACVR1B and WNT10A can be regulated by miR-133 and miR-23b, the protein levels of ACVR1B and WNT10A in sheep dermal fibroblasts transfected with miR-133 and miR-23b mimics and inhibitor were determined by Western blotting. The results showed that miR-133 and miR-23b mimics dramatically decreased the protein levels of ACVR1B and WNT10A, while miR-133 and miR-23b inhibitors increased the ACVR1B (Fig. 8f, Additional file 5: Fig S4) and WNT10A (Fig. 8g, Additional file 5: Fig S4) protein levels in sheep dermal fibroblasts. Furthermore, the results were in agreement with the decreased expression of ACVR1B and WNT10A (Fig. 8d, 8e).