Skip to main content

Integrated analysis of miRNAs and mRNA profiling reveals the potential roles of miRNAs in sheep hair follicle development

Abstract

Background

Merino sheep exhibit high wool production and excellent wool quality. The fleece of Merino sheep is predominantly composed of wool fibers grown from hair follicles (HFs). The HF is a complex biological system involved in a dynamic process governed by gene regulation, and gene expression is regulated by microRNAs (miRNAs). miRNA inhibits posttranscriptional gene expression by specifically binding to target messenger RNA (mRNA) and plays an important role in regulating gene expression, the cell cycle and biological development sequences. The purpose of this study was to examine mRNA and miRNA binding to identify key miRNAs and target genes related to HF development. This will provide new and important insights into fundamental mechanisms that regulate cellular activity and cell fate decisions within and outside of the skin.

Results

We analyzed miRNA data in skin tissues collected from 18 Merino sheep on four embryonic days (E65, E85, E105 and E135) and two postnatal days (D7 and D30) and identified 87 differentially expressed miRNAs (DE-miRNAs). These six stages were further divided into two longer developmental stages based on heatmap cluster analysis, and the results showed that DE-mRNAs in Stage A were closely related to HF morphogenesis. A coanalysis of Stage A DE-mRNAs and DE-miRNAs revealed that 9 DE-miRNAs and 17 DE-mRNAs presented targeting relationships in Stage A. We found that miR-23b and miR-133 could target and regulate ACVR1B and WNT10A. In dermal fibroblasts, the overexpression of miR-133 significantly reduced the mRNA and protein expression levels of ACVR1B. The overexpression of miR-23b significantly reduced the mRNA and protein expression levels of WNT10A.

Conclusion

This study provides a new reference for understanding the molecular basis of HF development and lays a foundation for further improving sheep HF breeding. miRNAs and target genes related to hair follicular development were found, which provided a theoretical basis for molecular breeding for the culture of fine-wool sheep.

Peer Review reports

Background

Subo Merino (SBM) sheep are breed reared in China that produce superfine wool. Wool fiber diameter is typically 17–19 μm, which is larger than the wool quality count of 80 [1] and greatly impacts the fine-wool sheep industry. SBM sheep are an essential element of animal husbandry in China. Hair follicles (HFs), which control sheep wool growth and development, are thin structures with a complex morphology and a periodic growth pattern [2]. In addition to the epithelial and hypodermal layers, hairs also include connective tissue sheaths, inner root sheaths, outer root sheaths, hair bulbs, and hair shafts [3]. The HFs of fine-wool sheep are composed of primary follicles (PFs) and secondary follicles (SFs); PFs appear early, whereas SFs appear late [4]. Compared with PFs, SFs have a more obvious effect on wool yield and quality. The larger the proportion of SFs, the smaller the wool fiber diameter is [5]. The most important factor affecting HF density is SF redifferentiation, which can be increased to improve wool fineness. It has been shown that HFs present distinct characteristics of gene expression and functions [6]. Merino wool follicles have been described in detail, and it is well established that no new follicles appear after birth [7,8,9,10,11,12]. Within 75 days of gestation, the first follicles begin to form in sheep fetus PFs, and they produce fibers by 90 days [7, 8, 13]. It is not until approximately 85 days of gestation that SFs begin to develop. At approximately 105 days of gestation, secondary-derived follicles (SDs) begin to branch off from SFs [14]. Because HFs fully mature before birth, their numbers do not increase after birth. A study by Zhao et al. [15] from our group investigated the morphogenesis of ovine HFs at six developmental stages based on histomorphological observations. They found that placodes and dermal condensate started to form at E65, indicating the induction of HFs. At E85, the number of PFs increased, and SFs began to form. SFs began to differentiate at E105, and the number of SDs began to increase at E135. HFs had matured and showed a complete structure at E135, and most of the hair shafts had emerged through the epidermis [15].

Hair bud formation in the epidermis is initiated by signals from the dermis, followed by the release of TGF-β [16], EDA/EDAR [17], NF-κB [18, 19], Noggin/Lef-1 [20], shh [21], BMP-2/4/7 [22], WNT/β-catenin [23,24,25,26,27], and FGF in hair buds, which induces dermal fibrogenic cells to form the dermal papillae that regulate keratinocyte activity [28]. Embryonic fibroblasts are essential for HF morphogenesis and wound closure, but they are subsequently unable to induce new HFs unless embryonic programs are reactivated. Fibroblasts are the main resident cells of the skin dermis and participate in embryonic HF morphogenesis and wound healing, since extracellular matrix is remodeled and collagen is deposited in the wound bed as a result of their activity [29]. Importantly, dermal fibroblasts originate from two distinct developmental lineages with unique functions that differentially mediate the response to epidermal signals such as Hedgehog signaling [30]. The cells of dermal papillae release a second signal to stimulate epithelial cell proliferation, which initiates differentiation and results in a complete HF structure. Considering the complexity of HFs, studies on fetal skin have been valuable for identifying all candidate genes that appear to be developmentally regulated.

The majority of genetic variants discovered in genome-wide association studies (GWAS) occur in noncoding RNAs [31]. Gene expression is regulated by noncoding RNAs, including miRNAs. MiRNAs, which are endogenously expressed in mammals [32,33,34], are 18–22-nucleotide (nt), small noncoding RNAs that are evolutionarily conserved and negatively regulate the expression of target mRNAs mainly by repressing translation or, in some cases, by cleaving mRNA transcripts [35, 36]. miRNAs play a key role in many biological processes by regulating gene expression at the posttranscriptional level. This will result in a change in the expression levels of the set of target genes by approximately 15–40% at the mRNA and protein levels [36]. miRNAs act as sponges to adsorb circular RNAs and serve as bridges for ceRNAs. miRNAs can also compete with lncRNAs for binding. They play a key role in gene expression regulatory networks. Individual miRNAs exhibit unique spatiotemporal expression patterns in developing and cycling HFs [36]. Therefore, it is important to understand the role of miRNAs in HF development. Several miRNAs have been identified from livestock species. However, compared with the results reported livestock such as goats (267 precursors, 436 mature [CHIR_2.0]) and cows (1064 precursors, 1025 mature [Btau_5.0.1]), the number of miRNAs identified in sheep is quite low (106 precursors, 153 mature [Oar_v4.0]), particularly in skin tissues. Recently, several studies have identified genes and miRNAs involved in HF development in sheep [37,38,39,40] and goats [41,42,43].

miR-205 and miR-203 are squamous epithelial miRNAs [44,45,46]. During skin development, Yi et al. showed that miR-205 is highly enriched in epithelial progenitors and stem cells. In addition, they found that miR-205 regulates PI(3)-kinase signaling and acts as an HF stem cell activator with a crucial role in HF morphogenesis [44]. Moreover, they showed that miRNA-203 is induced during stratification and differentiation in the skin [45]. miR-125b acts as a repressor of the differentiation of HF stem cells and is therefore crucial to the onset of anagen [47, 48]. The miRNAs miR-31 [49] and miR-214 [50] are abundantly expressed in proliferating hair matrix keratinocytes and the outer root sheath during anagen. During the hair cycle, miR-31 modulates gene expression programs in the skin and HFs. miR-31 modulates hair growth at least partly by modulating BMP and FGF signaling pathways [49]. miR-214 inhibits hair growth by inhibiting the Wnt and Shh signaling pathways [50]. miR-24 is involved in the regulation of differentiation programs during HF development and is predominantly expressed in differentiated keratinocytes of the inner root sheath [51]. When miR-24 is overexpressed in proliferating cells during skin morphogenesis, abnormal HF development is associated with reduced proliferation and premature differentiation [52]. HF and skin miRNA databases were enriched by these studies, and our understanding of miRNA regulation in skin and HFs was improved. However, the mechanisms by which miRNAs regulate HF development in fine-wool sheep are not well understood.

In this study, to elucidate the molecular mechanisms of HF development in sheep, we sequenced miRNAs across six important HF development stages (i.e., E65, E85, E105, E135, D7, D30) in the skin tissue of 18 Merino sheep. To understand the functions of miRNAs and their target genes in HFs, we constructed an miRNA–mRNA regulatory network including the results of previous mRNA analyses conducted by our research group. Additionally, the targeting relationships between miRNAs and mRNAs were verified with a dual-luciferase reporter gene system, and RT–qPCR and Western blotting verified the overexpression/inhibition of miRNAs in sheep skin dermal fibroblasts. This work provides new insight into the miRNA and mRNA interaction networks involved in HF development and will potentially benefit wool quality control in the wool industry.

Results

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 Ensembl sheep genome 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 the mapping of more than 80% of the clean reads in the 18 libraries (Additional file 1: Table S1), and the length distribution of small RNAs differed among the eighteen libraries. We filtered raw data, and the results showed that the coverage of clean miRNA reads (Q30) 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. S1b). The sequencing depth and coverage of the samples met the requirements of this experiment and indicated that they could be used for subsequent experiments.

RT–qPCR validation of miRNA-seq

We compared the RT–qPCR results with the miRNA-seq results (Fig. 1). We found by RT–qPCR that the skin expression levels of miR-148a in three periods (E135, D7 and D30) and miR-154a-5p in four periods (E65 to E135) were inconsistent with the sequencing results. The RT–qPCR 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.

Fig. 1
figure 1

Results of RT–qPCR and miRNA-seq. The x-axis represents the time period, and the y-axis represents relative expression

Identification of differentially expressed miRNAs

A total of 87 DE-miRNAs and 446 novel DE-miRNAs were identified in the six HF development stages in SBM sheep. Previous research by our group identified 7879 DE-mRNAs [53]. A total of 11 DE-miRNAs and 104 novel DE-miRNAs were identified in E65 vs. E85, among 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 identified. 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. The E65 vs. D30 comparison revealed the greatest number of DE-miRNAs (60), and E135 vs. D7 showed the fewest DE-miRNAs (only 9). These results 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, among which E65 vs. E85 showed 2 DE-miRNAs (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 miRNAs (miR-136, miR-380-5p, miR-411b-5p, etc.) were specifically expressed E105 vs. E135. Two miRNAs (miR-1197-5p and miR-410-5p) and one miRNA (miR-487a-3p) were specifically expressed in E135 vs. D7 and D7 vs. D30, respectively. E65vs. D30 showed the greatest number of specifically expressed miRNAs (25).

Fig. 2
figure 2

a The DE-miRNAs and novel DE-miRNAs between groups were obtained via pairwise comparisons; Up indicates a significantly upregulated miRNA; Down indicates a significantly downregulated miRNA. b Venn diagram of DE-miRNAs

Enrichment analysis of miRNA target genes

To further elucidate the functions of the DE-miRNAs, we also examined the top 10 GO functional enrichment biological process (BP) (Fig. 3), cellular component (CC) and molecular function (MF) terms (Additional file 3: Fig. S2) and KEGG pathways (Additional file 4: Fig. S3) of target genes in adjacent comparison groups. In E65 vs. E85, the enriched GO terms were mainly related to the positive regulation of vacuole organization, innate immune response, extrinsic apoptotic signaling pathway and negative regulation of the Notch signaling pathway. The GO terms of E85 vs. E105 were enriched in cell proliferation involved in outflow tract morphogenesis and segmentation. In E105 vs. E135, the terms were mainly concentrated in the regulation of the apoptotic process and the regulation of programmed cell death. In E135 vs. D7, the terms were mainly enriched in BPs related to the negative regulation of the Notch signaling pathway, positive regulation of vacuole organization, and outflow tract septum morphogenesis. In D7 vs. D30, the terms were mainly enriched in cell cycle checkpoint, tube closure and chordate embryonic development. 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 identified in E65 vs. E85 were directly related to cellular immunity and apoptosis. The target genes identified in E105 vs. E135 were significantly correlated with apoptosis. The target genes identified in E135 vs. D7 were related to kidney development, and those in E65 vs. D30 were related to cardiac and neurological development.

Fig. 3
figure 3

Top 10 GO classification statistics of miRNA target genes in different stages of hair follicle morphogenesis

We performed KEGG enrichment analysis and found that the target genes were mainly enriched in the proteasome, the AMPK signaling pathway, protein processing in the endoplasmic reticulum and other pathways. Similarly, they were enriched in the neomycin, kanamycin and gentamicin biosynthesis pathways (Additional file 4: Fig. S3).

DE-mRNA and DE-miRNA cluster analysis

To screen mRNAs and miRNAs according to heatmap analysis, the DE-mRNAs identified at the first two time points (E65 and E85) and the last four time points (E105, E135, D7 and D30) were combined as Stage A and Stage B DE-mRNAs, respectively (Fig. 4a). Similarly, the DE-miRNAs from the first three time points (E65, E85 and E105) and the last three time points (E135, D7 and D30) were combined as Stage A and Stage B DE-miRNAs, respectively (Fig. 4c). Then, we generated a Venn diagrams for Stage A and Stage B. The numbers of Stage A DE-mRNAs and DE-miRNAs were 1562 and 21, respectively, and the numbers of Stage B DE-mRNAs and DE-miRNAs were 2100 and 28, respectively. The total 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 the respective HF development stage (Stage A). The remaining 1672 DE-mRNAs are candidate genes associated with the HF 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).

Fig. 4
figure 4

a Heatmap of DE-mRNAs. b Venn diagram of DE-mRNAs. c Heatmap of DE-miRNAs. d Venn diagram of DE-miRNAs

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 the stage A DE-mRNAs (Fig. 5). In the BP category, the DE-mRNAs were significantly enriched in the negative regulation of the canonical Wnt signaling pathway, HF development, establishment of the skin barrier, skin development, the 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 BPs have been reported to be associated with HF morphogenesis, skin development or even hair formation. We visualized the genes involved in these BPs 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.

Fig. 5
figure 5

a Statistical map of GO classifications in Stage A. b Analysis of GO networks related to HF development (triangles represent GO terms, circles represent genes). c Statistical map of signaling pathways s in Stage A. d Analysis of KEGG networks related to hair follicle development (squares represent KEGG pathways, circles represent genes)

A list of significant pathways was obtained via the KEGG analysis of the 1134 DE-mRNAs in Stage A (Fig. 5c). Four significant pathways were identified: the Hippo signaling pathway, the PI3K-Akt signaling pathway, the TGF-β signaling pathway and the 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 identified 20 genes, including TGFβ2, WNT10A, PTCH1, BMP2, ACVR1B, WNT5A, WNT16, LAMA5, WNT10A, and WNT4, that may be involved in HF development (Fig. 5d). In summary, we identified 19 genes (LAMA5, ITGB4, ITGA3, THBS1, NOG, ACVR1B, INHBA, TGFβ2, BMP2, BMP5, WNT10A, SOX2, WNT16, WNT14, WNT5A, PTCH1, PTCH2, SHH, and FZD1) that play a key role in Stage A of HF development.

A list of significant pathways was obtained via the 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: the positive regulation of NF-kappa B transcription factor activity (Fig. 6a). Among the identified pathways, the PPAR signaling pathway has been best studied in the context of HF development (Fig. 6b). Our study suggested that the genes in this pathway are associated with the late stage of HF development and maturation.

Fig. 6
figure 6

a Statistical map of GO classifications in Stage B. b Diagram of signaling pathways in Stage B

Analysis of miRNA-mRNA and miRNA target DE-mRNA networks

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.

Fig. 7
figure 7

a DE-mRNA and DE-miRNA networks in Stage A. b Relative expression of miR-133 and ACVR1B in skin tissues at different stages of hair follicle morphogenesis. c Relative expression of miR-23b and WNT10A in skin tissues at different stages of hair follicle morphogenesis

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 maximal 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 evaluate the targeting relationships between these miRNAs and genes, we constructed two luciferase reporter vectors, psiCHECK2-WNT10A and psiCHECK2-ACVR1B, for the transfection of 293 T 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’-untranslated regions (UTRs) s (Fig. 8a). These reporters were cotransfected with miR-133 and miR-23b mimics or a negative control (mimic-NC) into 293 T cells. Luciferase reporter activity was measured in dual luciferase assays. As shown in Fig. 8, luciferase activity significantly decreased when cells were cotransfected with miR-133 and miR-23b mimics compared to the transfection of either mutant or empty controls. The overexpression of miR-133 and miR-23b mimics decreased the luciferase activity of psiCHECK2-ACVR1B (Fig. 8b) and psiCHECK2-WNT10A (Fig. 8c).

Fig. 8
figure 8

miR-133 could target and regulate ACVR1B, and miR-23b could target and regulate WNT10A. a The predicted sites of miR-133 binding to ACVR1B and miR-23b binding to WNT10A. b A luciferase reporter gene assay revealed that the overexpression of miR-133 significantly quenched wild-type ACVR1B fluorescence. c A luciferase reporter gene assay revealed that the overexpression of miR-23b significantly quenched wild-type WNT10A fluorescence. d RT–qPCR quantification indicated that the overexpression of miR-133 significantly reduced the mRNA expression level of ACVR1B. e RT–qPCR quantification indicated that the overexpression of miR-23b significantly reduced the mRNA expression level of WNT10A. f Western blot analysis indicated that the overexpression of miR-133 significantly reduced the protein expression level of ACVR1B. g Western blot analysis indicated that the overexpression of miR-23b significantly reduced the protein expression level of WNT10A. h Analysis of the relative protein expression of ACVR1B. i Analysis of the relative protein expression of WNT10A. **p < 0.01, *p < 0.05, ns p > 0.05

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 determine whether ACVR1B and WNT10A could 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 inhibitors were determined by RT–qPCR. The miR-133 and miR-23b mimics dramatically decreased the mRNA levels of ACVR1B and WNT10A, while miR-133 and miR-23b inhibitors increased 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 could 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 ACVR1B (Fig. 8f, h, Additional file 5: Fig S4) and WNT10A (Fig. 8g, i, 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, e).

Discussion

Microanatomy and cellular activity are dramatically altered during HF development and regeneration and are controlled by multiple signaling pathways, transcription factors and epigenetic regulators, including miRNAs. In the different HF cell lineages, miRNAs and their targets form remarkably diverse regulatory networks that play vital roles in gene expression programming. Andl and Botchkareva [36] summarized the roles of miRNAs in the control of HF development, cycling and hair pigmentation and presented future directions for this exciting and rapidly growing area of research. Through cell-specific or tissue-specific interference with miRNA biogenesis, it has been proven that miRNAs play indispensable roles in development and organogenesis. During embryonic development, the structural surface ectoderm-specific deletion of the miRNA processor Dicer or DgCr8 in the skin leads to serious abnormalities in HF development, characterized by the inability of HFs to invaginate into the dermis, leading to major structural defects and an inability to produce hair stems [3, 54]. Of course, miRNAs are essential not only for the development of HFs but also for postnatal HF growth. When Drosha or Dicer is deleted after HF initiation, fully developed HFs show a large number of apoptotic cells and structural changes in hair stems in the hair matrix, ultimately leading to the degradation of the HF [55]. These studies provide significant evidence that the global deletion of miRNAs in the embryonic epithelium significantly changes the molecular interaction between epithelial cells and dermal cells, resulting in the failure to form normal HFs, and show that miRNAs are necessary to control proliferation and apoptosis in adult HFs.

The skin sends signals to initiate hair bud formation, and the bud releases signaling factors that induce fibrotic cells to form dermal papillae that regulate keratinocyte activity [28]. Fibroblasts are major mesenchymal cells that deposit collagen and elastic fibers in the extracellular matrix (ECM); even within a single tissue, fibroblasts exhibit considerable functional diversity [29]. There are two distinct lines of fibroblasts that form the connective tissue of the skin. One forms the upper dermis, including the dermal papilla that regulates hair growth and the arrector pili muscle, which controls piloerection. The upper dermis is required for HF formation [28, 29].

The initiation of HFs involves a series of signaling pathways that link epidermal cells and the dermal papillae, such as the Wnt/β-catenin [56, 57], TGF-β [58] and Hippo signaling pathways. WNT, Hippo and TGF-β signals are required for the initiation of HF development. The current study also revealed that the above signaling molecules are expressed in sheep skin during the onset of secondary HF development. Based on the results of our study, we consider the DE-mRNAs and miRNAs of Stage A to have more practical significance in HF development than those of Stage B. E65, E85 and E105 represent the key periods for the formation and differentiation of PFs, SFs and SDs. We performed a combined analysis of the mRNA and miRNA sequencing results and ultimately identified 9 DE-miRNAs and 17 DE-mRNAs with targeting relationships in Stage A. Combined with existing research reports, our DE-mRNA enrichment analysis results revealed that the ACVR1B gene plays a key role in the TGF-β signaling pathway. WNT10A is expressed in the Wnt and Hippo signaling pathways. Therefore, we believe that ACVR1B and WNT10A play key roles in the initiation of HF development. Finally, we selected ACVR1B and WNT10A and their targeted miRNAs for subsequent validation. However, the dual luciferase reporter gene results showed that miR-23b and ACVR1B did not present a targeted regulatory relationship. The in vivo functions of ACVR1B have been difficult to study because ACVR1B knockout mice die during embryogenesis. Mice with conditional disruption of ACVR1B display various degrees of hairlessness at postnatal day 5, and the phenotype becomes more severe with age [29]. HFs that develop during morphogenesis are later disrupted by delays in hair cycle re-entry and the failure of cycling and regrowth of the hair shaft and the inner root sheath, resulting in severe hair loss [29]. This study demonstrates a specialized role of ACVR1B in hair cycling as well as HF development. Therefore, ACVR1B signaling is required for both HF development and cycling. In our study, we found that miR-133 could target and regulate ACVR1B. In dermal fibroblasts, the overexpression of miR-133 significantly reduced the mRNA and protein expression levels of ACVR1B.

Andl [59] et al. used a secreted inhibitor of essential WNT coreceptors and showed that signaling by canonical WNT proteins is required for the initiation of all types of HFs. These results demonstrate that the actions of canonical WNT proteins in the skin precede the localized expression of genes that regulate HF placode formation and indicate that canonical WNT signaling is required for the generation of, or the initial response to, the first dermal message. Wu et al. [60] showed that WNT10B can act via the Wnt/β-catenin signaling pathway to inhibit HF development in vitro. Although these Wnts are expressed in the epidermal component of HFs and are candidate mediators of this activity in vivo, it is possible that another Wnt acting through the same signal transduction pathways may normally serve this function during the hair cycle. Reddy et al. [61] found that FZD1 expression was correlated with WNT10A and WNT10B expression in the placode. They suggested that canonical WNT signaling in the placodes and the dermal condensate of developing HFs is likely activated by WNT10A and WNT10B expressed in and secreted from epithelial cells that then bind to FZD1 expressed in the placode epithelium and dermal condensate. WNT10A and WNT10B are continuously present in the follicular epithelium at later germ and bulbous peg stages. Similarly, they are strongly expressed in the epithelial cone surrounding the dermal papilla during the bulbous peg stage. However, Wnt signaling through the β-catenin pathway is sufficient to maintain, but not restore, the anagen-phase characteristics of dermal papilla cells [62]. It was shown that WNT10A is expressed in the outer root sheath, inner root sheath, matrix and hair shaft of anagen follicles in rats [57]. The characterization of the molecular pathways controlling differentiation and proliferation in mammalian HFs is central to our understanding of the regulation of normal hair growth, the foundations of hereditary hair loss diseases and the origins of follicle-based tumors. Millar et al. [63] demonstrated that the expression of WNT10A, WNT11, and β-catenin was increased at days 14 and 21 post wounding, when tissue remodeling occurred. In the epithelial tongue and neoepidermis, the WNT10A signal contributes to cutaneous wound repair by stimulating ECM maturation during the initial stage of healing [64]. WNT10A expression, especially in fibroblasts, can be crucial for achieving various potentially beneficial effects in wound healing [65]. Studies have also shown that WNT10A plays an important role in the pathogenesis of idiopathic pulmonary fibrosis via TGF-β activation and may be a sensitive predictor of the onset of an acute exacerbation of idiopathic pulmonary fibrosis [66]. In our study, we found that miR-23b could target and regulate WNT10A. In dermal fibroblasts, the overexpression of miR-23b significantly reduced the mRNA and protein expression levels of WNT10A.

Conclusion

In summary, we performed miRNA-seq of the miRNAs in six developmental stages of SBM HFs and conducted a bioinformatic analysis. We identified a total of 87 DE-miRNAs and 446 novel DE-miRNAs. We further performed mRNA–miRNA coanalysis in the six developmental stages of SBM HFs. We found that sheep HF development could be divided into two developmental stages (Stage A and Stage B). The mRNAs (E65 and E85) and miRNAs (E65 to E105) of stage A were closely related to HF morphogenesis. We found miRNAs and genes associated with hair follicular development, providing a theoretical basis for molecular breeding for the production of fine-wool sheep. Many more studies are needed to further delineate the significance of functional miRNA–mRNA interactions in HF biology. This research provides new insights into how the skin regulates cellular activity and cell fate decisions.

Methods

Animal selection and skin tissue preparation

SBM is a subtype of the Merino breed of sheep reared in China and is famous for the excellent quality and yield of its wool and its high survival rate. The tested individuals were selected from the sheep herd located at the Kechuang Animal Husbandry Breeding Center, Xinjiang, China. Twenty healthy SBM ewes were artificially inseminated with fresh sperm from the same SBM ram and then managed in the same flock. The day of insemination was designated embryonic day 0 (E0). Embryonic (E65, E85, E105 and E135) skin tissue collection has been described previously [53]. Skin tissue collection on postnatal days 7 and 30 (D7 and D30) has also been described previously [53]. Three biological replicates were generated for each of six developmental stages representing groups. All eighteen skin tissue samples were stored at − 80 °C.

Total RNA isolation and sequencing

Total RNA was isolated from the tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). For miRNA sequencing, the integrity of the total RNA samples was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., USA), and samples with RNA integrity numbers (RINs) above 7.0 were used for sequencing. Paired-end libraries were synthesized by using the QIAseq miRNA Library Kit (Qiagen, Germany) following the QIAseq miRNA Library Kit Guide. The products were then purified and enriched via PCR to create the final cDNA library. Purified libraries were quantified with a Qubit® 2.0 fluorometer (Life Technologies, USA). Clusters were generated by using cBot with the library diluted to 10 pM and were then sequenced on the Illumina HiSeq Xten platform (Illumina, USA).

Analysis of sequencing data

The sequence analysis of DE-mRNAs is detailed in another article we previously published [53]. miRNA-seq reads were filtered, and clean reads were mapped to the sheep (Ovis aries V3.1) mature miRNA database in miRBase (v21) and the Rfam database (ftp://selab.janelia.org/pub/Rfam) to match them with known long noncoding RNA, miRNA, ribosomal RNA, small nuclear RNA, small nucleolar RNA, and transfer RNA sequences. miRNA abundance was expressed as counts of exon model per million mapped reads (CPM). DE-miRNAs with |log2 (FC)| value > 1 and p value < 0.05 were considered significantly modulated. Clean reads that were not mapped to sheep mature miRNAs were then mapped to the Oar 3.1 version of the sheep genome sequence using Bowtie software, and the distribution of mapped reads on the chromosomes was calculated. Novel miRNAs were also predicted by using miRDeep2. In the mixed sequencing data of all samples of the studied species, various non-miRNA sequences were filtered out, and the remaining sequencing data (defined as clean miRNA reads) were used for the identification and structural prediction of new miRNAs. Predicting new miRNAs requires reference genome sequences. In this program, miRDeep software is used in combination with homologous miRNA sequences of closely related species, including the application of RNA secondary structure prediction software such as RNAfold, to predict and identify the mature forms (Star miRNA and mature miRNA) and precursors of new miRNAs. For each newly predicted miRNA, miRDeep software applies various evaluation methods to evaluate the prediction accuracy. The expression patterns of DE-miRNAs were measured by systematic cluster analysis to explore the similarities and compare the relationships between the different libraries.

RT–qPCR validation of miRNA-seq results

To confirm the differential expression of genes revealed by miRNA-seq, a total of nine miRNAs identified as differentially expressed among different developmental stages were randomly chosen for RT–qPCR validation. Total RNA was extracted using TRIzol reagent (Invitrogen). A Poly(A) Tailing Kit (Ambion, Waltham, MA, USA) was used to add poly(A) tails to miRNAs. We used the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) and gene-specific primers or random primers to generate cDNA. RT–qPCR was performed in a CFX96™ Real-Time System (Bio-Rad, USA) using the miRcute Plus miRNA qPCR Kit (SYBR Green) (Tiangen, China). U6 snRNA was employed as an endogenous control for miRNA. In the miRNA RT‒qPCR assay, the thermal cycling conditions were 95 °C for 15 min, followed by 45 cycles at 94 °C for 20 s and 60 °C for 34 s. Three biological and technical replicates were conducted. The 2−ΔΔCt method was used to determine relative miRNA abundance. Primer sequences are displayed in Additional file 6: Table S2.

GO enrichment and KEGG pathway analysis

The GO (http://geneontology.org) and KEGG (http://www.genome.ad.jp/kegg/) databases of DAVID were used to perform functional annotations. DE-mRNAs identified in Stages A and B of HF development were analyzed. Based on the false discovery rate (FDR), we determined the GO terms and metabolic pathways significantly associated with the gene lists. An FDR ≤ 0.05 was applied to significant candidate genes associated with BPs or metabolic pathways.

miRNA target gene prediction and construction of an interaction network

In animals, miRNA binds to the 3’- UTR of the target gene mRNA in an incomplete complementary pairing manner. The genes in Stage A related to HF morphogenesis obtained from the above enrichment analysis were predicted using the miRanda, TargetScan (http://www.targetscan.org), PicTar (https://pictar.mdc-berlin.de), RNAhybrid (https://bibiserv.cebitec.uni-biele), miRWalk (http://mirwalk.umm.uni-heidelberg.de/), RNA22 (http://cbcsrv.watson.ibm.com/rna22.html) and starBase (http://starbase.sysu.edu.cn/) tools to predict miRNA target genes. Then, more than three cross-target genes were selected as target mRNAs in the Venn analysis. mRNAs whose miRNA-targeted mRNAs overlapped with DE-mRNAs in phase A were selected as DE-mRNAs. Finally, Cytoscape was used for the visualization and analysis of the DE-mRNA and DE-miRNA networks.

Cell culture and transfection

HEK-293 T cells were used for the dual-luciferase reporter gene assay and cultured at 37 °C in Dulbecco’s modified Eagle’s medium (DMEM) (Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS) (Invitrogen, Carlsbad, CA, USA), 1.5 mM l-glutamine (Invitrogen, Carlsbad, CA, USA), 100 U/mL penicillin (Invitrogen, Carlsbad, CA, USA), and 100 mg/mL streptomycin (Invitrogen, Carlsbad, CA, USA) in a humidified incubator in an atmosphere containing 5% CO2 (Thermo, Waltham, MA, USA). We seeded the cells into 24-well plates, and when the cells reached 80% confluence, we used Lipofectamine™ 3000 (Invitrogen, CA) to cotransfect 293 T cells with the constructed double luciferase reporter gene vector (psiCHECK2-WNT10A-WT, psiCHECK2-WNT10A-MUT, psiCHECK2-ACVR1B-WT and psiCHECK2-ACVR1B-MUT) along with miR-23b and miR-133 (mimic and mimic-NC).

Sheep dermal fibroblasts were used to verify the effects of miR-23b and miR-133 on their targeted mRNAs. The skin tissue of healthy postnatal lambs at 7 days of age was collected at the Xinjiang Kechuang Breeding Center. Sheep skin was treated and cultured at 37 °C with 5% CO2, and the medium was then changed every three days to observe cell growth in the culture dish to obtain primary fibroblasts. Cells were passaged when the cell confluency was more than 90%. When the cells reached to the fourth generation, they were inoculated into a 6-well plate and grown until the confluence rate reached 80%. We used Lipofectamine™ 3000 (Invitrogen, CA)-transfected fibroblasts with miRNA overexpression (miR-23b mimic and miR-133 mimic) or knockout (miR-23b inhibitor and miR-133 inhibitor) or a negative control (mimic-NC and inhibitor-NC) according to the manufacturer's instructions.

Dual luciferase reporter gene assay

Luciferase reporters used to verify the targeting relationship between miRNAs and mRNAs were generated based on the psiCHECK2 vector (Promega). To construct psiCHECK2-WNT10A and psiCHECK2-ACVR1B, the complete 3’-UTRs of sheep WNT10A mRNA (676 nt, GenBank accession no. XM_004004935) and ACVR1B mRNA (2966 nt, GenBank accession no. XM_027967297), containing the putative miR-133 and miR-23b binding sites, were amplified and cloned into the psiCHECK2 vector. The luciferase reporter was cotransfected with miR-133 mimics, miR-23b mimics, and mimics-NC into 293 T cells by using Lipofectamine™ 3000 according to the manufacturer’s guidelines. Relative luciferase activity was measured with a Dual-Luciferase Reporter Assay System (Promega) and an Infinate M200 PRO microplate reader (Tecan, Shanghai, China).

RT–qPCR and Western blotting of fibroblasts

RNA was extracted using an RNAsimple Total RNA Kit (Tiangen, China) and reverse transcribed into cDNA using a FastQuant RT Kit (Tiangen, China) after sheep dermal fibroblasts had been transfected for 36 h. RT–qPCR was used to quantify WNT10A and ACVR1B mRNA expression in the different treatment groups using TB Green Premix Ex Taq II (TaKaRa, Dalian, China) following the manufacturer's instructions. The RT‒qPCR experiments were conducted in triplicate, and the relative abundance of mRNA was determined using the 2−ΔΔCt method. ANOVA was used for statistical comparison of more than two groups. All statistical analyses were conducted using SPSS 20.0 software (IBM Corp, Armonk, NY, USA). P < 0.05 was considered statistically significant.

For Western blot analysis, dermal fibroblast cells were lysed using RIPA Lysis Buffer II (Sangon, Shanghai, China). Protein concentrations in different groups were measured with a Pierce™ Rapid Gold BCA Protein Assay Kit (Thermo Scientific, Waltham, MA, USA). The proteins were separated by SDS–PAGE, transferred to polyvinylidene difluoride membranes, and probed with 1:1000 rabbit anti-WNT10A (ABclonal, A19090, China), 1:1000 rabbit anti-ACVR1B (ABclonal, A2279, China) and 1:1000 rabbit anti-GAPDH (ABclonal, A19056, China) antibodies and then with 1:1000 HRP goat anti-rabbit IgG (H + L) (ABclonal, AS014, China)-conjugated antibodies. All steps were carried out in accordance with the provided instructions. Signals were detected with an enhanced chemiluminescence (ECL) substrate kit (UltraSensitive) (Biosharp, Beijing, China). Proteins were detected and analyzed with the Wes Automated Western Blot Analysis System (Protein Simple, CA, US).

Availability of data and materials

All miRNA-seq data generated in this study were submitted to the NCBI SRA database under BioProject No. PRJNA705552 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA705552). All RNA-seq data generated in this study were submitted to the NCBI SRA database under BioProject No. PRJNA705554 ( https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA705554).

Abbreviations

HF:

Hair follicle

miRNA:

MicroRNA

mRNA:

Messenger RNA

DE-miRNA:

Differentially expressed miRNA

DE-mRNA:

Differentially expressed mRNA

SBM:

Subo merino

PF:

Primary follicle

SF:

Secondary follicle

SD:

Secondary-derived follicle

TGF-β:

Transforming growth factor beta

BMP:

Bone morphogenetic protein

RT–qPCR:

Reverse transcription quantitative real-time PCR

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

BP:

Biological processes

ECM:

Extracellular matrix

CPM:

Counts of exon model per million mapped reads

FC:

Fold-change

DAVID:

Database for Annotation, Visualization, and Integrated Discovery

FDR:

False discovery rate

3’ UTR:

3’-Untranslated region

References

  1. Sulayman A, Tian K, Huang X, Tian Y, Xu X, Fu X, et al. Genome-wide identification and characterization of long non-coding RNAs expressed during sheep fetal and postnatal hair follicle development. Sci Rep. 2019;9(1):8501. https://doi.org/10.1038/s41598-019-44600-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Galbraith H. Fundamental hair follicle biology and fine fibre production in animals. Animal. 2010;4(9):1490–509. https://doi.org/10.1017/S175173111000025X.

    Article  CAS  PubMed  Google Scholar 

  3. Andl T, Murchison EP, Liu F, Zhang Y, Yunta-Gonzalez M, Tobias JW, et al. The miRNA-processing enzyme dicer is essential for the morphogenesis and maintenance of hair follicles. Curr Biol. 2006;16(10):1041–9. https://doi.org/10.1016/j.cub.2006.04.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Wang J, Sui J, Mao C, Li X, Chen X, Liang C, et al. Identification of key pathways and genes related to the development of hair follicle cycle in cashmere goats. Genes-Basel. 2021;12(2):180. https://doi.org/10.3390/genes12020180.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Zhao B, Fu X, Tian K, Huang X, Di J, Bai Y, et al. Identification of SNPs and expression patterns of FZD3 gene and its effect on wool traits in Chinese merino sheep (Xinjiang Type). J Integr Agr. 2019;18(10):2351–60. https://doi.org/10.1016/S2095-3119(19)62735-8.

    Article  CAS  Google Scholar 

  6. Wang Q, Qu J, Li Y, Ji D, Zhang H, Yin X, et al. Hair follicle stem cells isolated from newborn yangtze river delta white goats. Gene. 2019;698:19–26. https://doi.org/10.1016/j.gene.2019.02.052.

    Article  CAS  PubMed  Google Scholar 

  7. Hardy H, Lyne A. The pre-natal development of wool follicles in merino sheep. Aust J Biol Sci. 1956;9(3):423–41.

    Article  Google Scholar 

  8. Brook AH, Short BF, Lyne AG. Formation of new wool follicles in the adult sheep. Nature. 1960;185:51. https://doi.org/10.1038/185051a0.

    Article  CAS  PubMed  Google Scholar 

  9. Chapman RE, Hopkins PS, Thorburn GD. The Effects of Fetal Thyroidectomy and Thyroxine Administration On the Development of Skin and Wool Follicles of Sheep Fetuses. J Anat. 1974; 117(Pt 2):419–32. https://www.ncbi.nlm.nih.gov /pmc/articles/PMC1231415/pdf/janat00385–0188.

  10. Hutchison G, Mellor DJ. Effects of maternal nutrition on the initiation of secondary wool follicles in foetal sheep. J Comp Pathol. 1983;93(4):577–83. https://doi.org/10.1016/0021-9975(83)90064-6.

    Article  CAS  PubMed  Google Scholar 

  11. Chapman RE, Hardy MH. Effects of intradermally injected and topically applied mouse epidermal growth factor on wool growth, skin and wool follicles of merino sheep. Aust J Biol Sci. 1988;41(2):261–8. https://doi.org/10.1071/bi9880261.

    Article  CAS  PubMed  Google Scholar 

  12. Rahmani W, Abbasi S, Hagner A, Raharjo E, Kumar R, Hotta A, et al. Hair follicle dermal stem cells regenerate the dermal sheath, repopulate the dermal papilla, and modulate hair type. Dev Cell. 2014;31(5):543–58. https://doi.org/10.1016/j.devcel.2014.10.022.

    Article  CAS  PubMed  Google Scholar 

  13. Priestley GC, Ryder ML, Head KW. Regeneration of wool follicles in autografts of sheep skin. Res Vet Sci. 1977;23(3):303–9. https://doi.org/10.1016/s0034-5288(18)33122-9.

    Article  CAS  PubMed  Google Scholar 

  14. Hocking EJ. Reduction in wool follicles prior to birth in merino sheep. Reprod Fert Develop. 1999;11(4–5):229. https://doi.org/10.1071/rd99049.

    Article  Google Scholar 

  15. Zhao B, Luo H, He J, Huang X, Chen S, Fu X, et al. Comprehensive transcriptome and methylome analysis delineates the biological basis of hair follicle development and wool-related traits in merino sheep. BMC Biol. 2021;19(1):197. https://doi.org/10.1186/s12915-021-01127-9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Oshimori N, Fuchs E. Paracrine TGF-beta signaling counterbalances BMP-mediated repression in hair follicle stem cell activation. Cell Stem Cell. 2012;10(1):63–75. https://doi.org/10.1016/j.stem.2011.11.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Wu Z, Wang Y, Han W, Yang K, Hai E, Ma R, et al. EDA and EDAR expression at different stages of hair follicle development in cashmere goats and effects on expression of related genes. Arch Anim Breed. 2020;63(2):461–70. https://doi.org/10.5194/aab-63-461-2020.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Zhang Y, Tomann P, Andl T, Gallant NM, Huelsken J, Jerchow B, et al. Reciprocal requirements for EDA/EDAR/NF-κB and Wnt/β-catenin signaling pathways in hair follicle induction. Dev Cell. 2009;17(1):49–61 (https://pubmed.ncbi.nlm.nih.gov/343207/).

    Article  Google Scholar 

  19. Laurikkala J, Pispa J, Jung H, Nieminen P, Mikkola M, Wang X, et al. Regulation of hair follicle development by the TNF signal ectodysplasin and its receptor edar. Development (Cambridge). 2002;129(10):2541–53. https://doi.org/10.1242/dev.129.10.2541.

    Article  CAS  Google Scholar 

  20. Botchkarev VA, Botchkareva NV, Nakamura M, Huber O, Funa K, Lauster R, et al. Noggin is required for induction of the hair follicle growth phase in postnatal skin. FASEB J. 2001;15(12):2205–14. https://doi.org/10.1096/fj.01-0207com.

    Article  CAS  PubMed  Google Scholar 

  21. Woo W, Zhen HH, Oro AE. Shh maintains dermal papilla identity and hair morphogenesis via a noggin-shh regulatory loop. Gene Dev. 2012;26(11):1235–46. https://doi.org/10.1101/gad.187401.112.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Rendl M, Polak L, Fuchs E. BMP signaling in dermal papilla cells is required for their hair follicle-inductive properties. Gene Dev. 2008;22(4):543–57. https://doi.org/10.1101/gad.1614408.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Lee YJ, Park SH, Park HR, Lee Y, Kang H, Kim JE. Mesenchymal stem cells antagonize IFN-induced proinflammatory changes and growth inhibition effects via Wnt/beta-Catenin and JAK/STAT pathway in human outer root sheath cells and hair follicles. Int J Mol Sci. 2021;22(9):4581. https://doi.org/10.3390/ijms22094581.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Du KT, Deng JQ, He XG, Liu ZP, Peng C, Zhang MS. MiR-214 regulates the human hair follicle stem cell proliferation and differentiation by targeting EZH2 and Wnt/beta-catenin signaling way in vitro. Tissue Eng Regen Med. 2018;15(3):341–50. https://doi.org/10.1007/s13770-018-0118-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Zhao B, Chen Y, Yang N, Chen Q, Bao Z, Liu M, et al. MiR-218-5P regulates skin and hair follicle development through Wnt/β-catenin signaling pathway by targeting SFRP2. J Cell Physiol. 2019;234(11):20329–41. https://doi.org/10.1002/jcp.28633.

    Article  CAS  PubMed  Google Scholar 

  26. Zhou L, Wang H, Jing J, Yu L, Wu X, Lu Z. Morroniside regulates hair growth and cycle transition via activation of the Wnt/β-catenin signaling pathway. Sci Rep-Uk. 2018;8(1):1–13. https://doi.org/10.1038/s41598-018-32138-2.

    Article  CAS  Google Scholar 

  27. Jin H, Zou Z, Chang H, Shen Q, Liu L, Xing D. Photobiomodulation therapy for hair regeneration: a synergetic activation of beta-CATENIN in hair follicle stem cells by ROS and paracrine WNTs. Stem Cell Rep. 2021;16(6):1568–83. https://doi.org/10.1016/j.stemcr.2021.04.015.

    Article  CAS  Google Scholar 

  28. Telerman SB, Rognoni E, Sequeira I, Pisco AO, Lichtenberger BM, Culley OJ, et al. Dermal Blimp1 acts downstream of epidermal TGFbeta and Wnt/beta-Catenin to regulate hair follicle formation and growth. J Invest Dermatol. 2017;137(11):2270–81. https://doi.org/10.1016/j.jid.2017.06.015.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Driskell RR, Lichtenberger BM, Hoste E, Kretzschmar K, Simons BD, Charalambous M, et al. Distinct fibroblast lineages determine dermal architecture in skin development and repair. Nature. 2013;504(7479):277–81. https://doi.org/10.1038/nature12783.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Frech S, Forsthuber A, Korosec A, Lipp K, Kozumov V, Lichtenberger BM. Hedgehog signaling in papillary fibroblasts is essential for hair follicle regeneration during wound healing. J Invest Dermatol. 2021. https://doi.org/10.1016/j.jid.2021.11.026.

    Article  PubMed  Google Scholar 

  31. Hindorff LA, Sethupathy P, Junkins HA, Ramos EM, Mehta JP, Collins FS, et al. Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proc Natl Acad Sci U S A. 2009;106(23):9362–7. https://doi.org/10.1073/pnas.0903103106.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Cai B, Zheng Y, Ma S, Xing Q, Wang X, Yang B, et al. Long non-coding rna regulates hair follicle stem cell proliferation and differentiation through PI3K/AKT signal pathway. Mol Med Rep. 2018;17(4):5477–83. https://doi.org/10.3892/mmr.2018.8546.

    Article  CAS  PubMed  Google Scholar 

  33. Lei Z, Sun W, Guo T, Li J, Zhu S, Lu Z, et al. Genome-wide selective signatures reveal candidate genes associated with hair follicle development and wool shedding in sheep. Genes-Basel. 2021;12(12):1924. https://doi.org/10.3390/genes12121924.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhang X, Bao Q, Jia C, Li C, Chang Y, Wu X, et al. Genome-wide detection and sequence conservation analysis of long non-coding RNA during hair follicle cycle of Yak. BMC Genomics. 2020;21(1):681. https://doi.org/10.1186/s12864-020-07082-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Han X, Chang L, Qiu Z, Lin M, Wang Y, Liu D, et al. Micro-injury induces hair regeneration and vitiligo repigmentation through Wnt/Β-Catenin pathway. Stem Cells Dev. 2022. https://doi.org/10.1089/scd.2021.0276.

    Article  PubMed  Google Scholar 

  36. Andl T, Botchkareva NV. MicroRNAs (miRNAs) in the control of HF development and cycling: the next frontiers in hair research. Exp Dermatol. 2015;24(11):821–6. https://doi.org/10.1111/exd.12785.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Yang H, Yang H, Shi G, Shen M, Yang JQ, Yang Y, et al. Expression profile analysis of microRNAs during hair follicle development in the sheep foetus. Biosci Biotechnol Biochem. 2019;83(6):1045–61. https://doi.org/10.1080/09168451.2019.1591261.

    Article  CAS  PubMed  Google Scholar 

  38. Lv X, Sun W, Yin J, Ni R, Su R, Wang Q, et al. An Integrated analysis of MicroRNA and mRNA expression profiles to identify RNA expression signatures in lambskin hair follicles in Hu Sheep. PLoS ONE. 2016;11(7):e157463. https://doi.org/10.1371/journal.pone.0157463.

    Article  CAS  Google Scholar 

  39. Lv X, Chen W, Wang S, Cao X, Yuan Z, Getachew T, et al. Integrated hair follicle profiles of microRNAs and mRNAs to reveal the pattern formation of hu sheep lambskin. Genes-Basel. 2022;13(2):342. https://doi.org/10.3390/genes13020342.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Gao W, Sun W, Yin J, Lv X, Bao J, Yu J, et al. Screening candidate microRNAs (miRNAs) in different lambskin hair follicles in hu sheep. PLoS ONE. 2017;12(5):e176532. https://doi.org/10.1371/journal.pone.0176532.

    Article  CAS  Google Scholar 

  41. Ge W, Zhang W, Zhang Y, Zheng Y, Li F, Wang S, et al. A single-cell transcriptome atlas of cashmere goat hair follicle morphogenesis. Genom Proteom Bioinform. 2021;19(3):437–51. https://doi.org/10.1016/j.gpb.2021.07.003.

    Article  Google Scholar 

  42. Hai E, Han W, Wu Z, Ma R, Shang F, Wang M, et al. Chi-miR-370–3p Regulates Hair Follicle Development of Inner Mongolian Cashmere Goats. G3 (Bethesda, Md.). 2021. https://doi.org/10.1093/g3journal/jkab091.

  43. Gao Y, Wang X, Yan H, Zeng J, Ma S, Niu Y, et al. Comparative transcriptome analysis of fetal skin reveals key genes related to hair follicle morphogenesis in cashmere goats. PLoS ONE. 2016;11(3):e151118. https://doi.org/10.1371/journal.pone.0151118.

    Article  CAS  Google Scholar 

  44. Yi R, O’Carroll D, Pasolli HA, Zhang Z, Dietrich FS, Tarakhovsky A, et al. Morphogenesis in skin is governed by discrete sets of differentially expressed microRNAs. Nat Genet. 2006;38(3):356–62. https://doi.org/10.1038/ng1744.

    Article  CAS  PubMed  Google Scholar 

  45. Yi R, Poy MN, Stoffel M, Fuchs E. A skin microRNA promotes differentiation by repressing ‘stemness.’ Nature. 2008;452(7184):225–9. https://doi.org/10.1038/nature06642.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Yu J, Peng H, Ruan Q, Fatima A, Getsios S, Lavker RM. MicroRNA-205 promotes keratinocyte migrationvia the lipid phosphatase SHIP2. FASEB J. 2010;24(10):3950–9. https://doi.org/10.1096/fj.10-157404.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Zhang L, Stokes N, Polak L, Fuchs E. Specific microRNAs are preferentially expressed by skin stem cells to balance self-renewal and early lineage commitment. Cell Stem Cell. 2011;8(3):294–308. https://doi.org/10.1016/j.stem.2011.01.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Qu H, Wu S, Li J, Ma T, Li J, Xiang B, et al. MiR-125b Regulates the Differentiation of Hair Follicles in Fine-wool Sheep and Cashmere Goats by Targeting MXD4 and FGFR2. Anim Biotechnol. 2021:1–08. https://doi.org/10.1080/10495398.2021.1968884.

  49. Mardaryev AN, Ahmed MI, Vlahov NV, Fessing MY, Gill JH, Sharov AA, et al. Micro-RNA-31 controls hair cycle-associated changes in gene expression programs of the skin and hair follicle. Faseb J. 2010;24(10):3869–81. https://doi.org/10.1096/fj.10-160663.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Ahmed MI, Alam M, Emelianov VU, Poterlowicz K, Patel A, Sharov AA, et al. MicroRNA-214 controls skin and hair follicle development by modulating the activity of the wnt pathway. J Cell Biol. 2014;207(4):549–67. https://doi.org/10.1083/jcb.201404001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Amelio I, Lena AM, Bonanno E, Melino G, Candi E. MiR-24 affects hair follicle morphogenesis targeting Tcf-3. Cell Death Dis. 2013;4(11):e922. https://doi.org/10.1038/cddis.2013.426.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Liu F, Zhang X, Peng Y, Zhang L, Yu Y, Hua P, et al. MiR-24 controls the regenerative competence of hair follicle progenitors by targeting Plk3. Cell Rep. 2021;35(10):109225. https://doi.org/10.1016/j.celrep.2021.109225.

    Article  CAS  PubMed  Google Scholar 

  53. He J, Zhao B, Huang X, Fu X, Liu G, Tian Y, et al. Gene network analysis reveals candidate genes related with the hair follicle development in sheep. BMC Genomics. 2022;23(1):428. https://doi.org/10.1186/s12864-022-08552-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Yi R, Pasolli HA, Landthaler M, Hafner M, Ojo T, Sheridan R, et al. DGCR8-dependent microRNA biogenesis is essential for skin development. Proc Natl Acad Sci. 2009;106(2):498–502. https://doi.org/10.1073/pnas.0810766105.

    Article  PubMed  Google Scholar 

  55. Teta M, Choi YS, Okegbe T, Wong G, Tam OH, Chong MMW, et al. Inducible deletion of epidermalDicer and drosha reveals multiple functions for miRNAs in postnatal skin. Development. 2012;139(8):1405–16. https://doi.org/10.1242/dev.070920.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Zhang R, Li Y, Jia K, Xu X, Li Y, Zhao Y, et al. Crosstalk between androgen and Wnt/β-catenin leads to changes of wool density in FGF5-knockout sheep. Cell Death Dis. 2020;11(5):407. https://doi.org/10.1038/s41419-020-2622-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Lin C, Yuan Y, Chen X, Li H, Cai B, Liu Y, et al. Expression of Wnt/β-catenin signaling, stem-cell markers and proliferating cell markers in rat whisker hair follicles. J Mol Histol. 2015;46(3):233–40. https://doi.org/10.1007/s10735-015-9616-5.

    Article  CAS  PubMed  Google Scholar 

  58. Massagué J. How cells read TGF-β signals. Nat Rev Mol Cell Bio. 2000;1(3):169–78. https://doi.org/10.1038/35043051.

    Article  Google Scholar 

  59. Andl T, Reddy ST, Gaddapara T, Millar SE. WNT signals are required for the initiation of hair follicle development. Dev Cell. 2002;2(5):643–53. https://doi.org/10.1016/S1534-5807(02)00167-3.

    Article  CAS  PubMed  Google Scholar 

  60. Wu Z, Zhu Y, Liu H, Liu G, Li F. Wnt10b promotes hair follicles growth and dermal papilla cells proliferation via Wnt/beta-catenin signaling pathway in rex rabbits. Biosci Rep. 2020;40(2):BSR20191248. https://doi.org/10.1042/BSR20191248.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Reddy ST, Andl T, Lu MM, Morrisey EE, Millar SE. Expression of frizzled genes in developing and postnatal hair follicles. J Invest Dermatol. 2004;123(2):275–82. https://doi.org/10.1111/j.0022-202X.2004.23215.x.

    Article  CAS  PubMed  Google Scholar 

  62. Shimizu H, Morgan BA. Wnt signaling through the β-catenin pathway is sufficient to maintain, but not restore, anagen-phase characteristics of dermal papilla cells. J Invest Dermatol. 2004;122(2):239–45. https://doi.org/10.1046/j.0022-202X.2004.22224.x.

    Article  CAS  PubMed  Google Scholar 

  63. Millar SE, Willert K, Salinas PC, Roelink H, Nusse R, Sussman DJ, et al. WNT signaling in the control of hair growth and structure. Dev Biol. 1999;207(1):133–49. https://doi.org/10.1006/dbio.1998.9140.

    Article  CAS  PubMed  Google Scholar 

  64. Bukowska J, Walendzik K, Kopcewicz M, Cierniak P, Gawronska-Kozak B. Wnt signaling and the transcription factor Foxn1 contribute to cutaneous wound repair in Mice. Connect Tissue Res. 2021;62(2):238–48. https://doi.org/10.1080/03008207.2019.1688314.

    Article  CAS  PubMed  Google Scholar 

  65. Wang K, Yamada S, Izumi H, Tsukamoto M, Nakashima T, Tasaki T, et al. Critical in vivo roles of WNT10A in wound healing by regulating collagen expression/synthesis in WNT10A-deficient Mice. PLoS ONE. 2018;13(3):e195156. https://doi.org/10.1371/journal.pone.0195156.

    Article  CAS  Google Scholar 

  66. Oda K, Yatera K, Izumi H, Ishimoto H, Yamada S, Nakao H, et al. Profibrotic Role of WNT10A Via TGF-beta signaling in idiopathic pulmonary fibrosis. Respir Res. 2016;17(1):39. https://doi.org/10.1186/s12931-016-0357-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the constructive comments from three reviewers.

Funding

This work was supported by the China Agriculture Research System (Nos. CARS-39) and the Agricultural Scientific and Technological Innovation Project of Shandong Academy of Agricultural Sciences (Nos. CXGC2022E10). The funders played no role in the study design, collection, analysis, data interpretation, manuscript writing, or decision to submit the manuscript for publication.

Author information

Authors and Affiliations

Authors

Contributions

KT conceived and designed the study. JH, XH, BZ and GL collected samples and/or generated data. JH, JM, and CW conducted the statistical analyses. GZ and YT managed the experimental sheep populations. JH, BZ and KT drafted and revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Kechuan Tian.

Ethics declarations

Ethics approval and consent to participate

The study was conducted in accordance with the Basel Declaration, and the report follows ARRIVE guidelines. All experimental protocols were approved by the Ethics Committee of the Institute of Animal Science and Veterinary Medicine of Shandong Academy of Agricultural Sciences (approval number: SAAS-2022-G24), and all tissue samples were collected according to the guidelines of the Ethics Committee. All embryo samples were obtained by cesarean section after 12 ewes were euthanized (exsanguination after electric shock syncope). Skin tissues were collected immediately after euthanasia and used for scientific research. On postnatal days D7 and D30, lambs were observed to lose consciousness, and their skin was collected in vivo from postnatal lambs, with a depth of approximately 2 cm2 × 3 mm. Afterward, the skin wounds were sprinkled with anti-inflammatory powder and carefully sutured. Until the wound had healed, we treated the lambs with particular care. We made every effort to minimize animal suffering and discomfort.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Result Statistics of miRNA Genome Comparison.

Additional file 2: Figure S1.

The distribution of the read from the sequencing data. (a) Read redundancy statistics of different samples. (b) Read length distribution of different samples.

Additional file 3: Figure S2.

Analysis of miRNA target genes go enrichment at different stages of hair follicle development. (a) GO enrichment analysis of E65 vs. E85 target genes. (b) GO enrichment analysis of E85 vs. E105 target genes. (c) GO enrichment analysis of E105 vs. E135 target genes. (d) GO enrichment analysis of E135 vs. D7 target genes. (e) GO enrichment analysis of D7 vs. D30 target genes. (f) GO enrichment analysis of E65 vs. D30 target genes. where red is BP (biological processes), green is CC (Cellular components) and blue is MF (Molecular functions).

Additional file 4: Figure S3.

Enrichment analysis of miRNA target gene KEGG pathway at different stages of hair follicle development.

Additional file 5: Figure S4.

(a) Western blot level of ACVR1B (target strip in red box dimension). (b) Western blot level of GAPDH (target strip in red box dimension). (c) Western blot level of WNT10A (target strip in red box dimension). (d) Western blot level of ACVR1B (target strip in red box dimension). All the images of western blots are cut prior to hybridization with antibodies. Because the markers are marked on both sides of the glue, the notch in the upper right corner (upper left corner) is to mark the front and back.

Additional file 6: Table S2.

Primer sequence.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

He, J., Huang, X., Zhao, B. et al. Integrated analysis of miRNAs and mRNA profiling reveals the potential roles of miRNAs in sheep hair follicle development. BMC Genomics 23, 722 (2022). https://doi.org/10.1186/s12864-022-08954-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08954-2

Keywords

  • Hair follicle
  • Merino sheep
  • miRNA–mRNA
  • Dermal fibroblasts
  • Western blot