Skip to main content

Comprehensive analysis of the expression profiles of mRNA, lncRNA, circRNA, and miRNA in primary hair follicles of coarse sheep fetal skin



The Qinghai Tibetan sheep, a local breed renowned for its long hair, has experienced significant deterioration in wool characteristics due to the absence of systematic breeding practices. Therefore, it is imperative to investigate the molecular mechanisms underlying follicle development in order to genetically enhance wool-related traits and safeguard the sustainable utilization of valuable germplasm resources. However, our understanding of the regulatory roles played by coding and non-coding RNAs in hair follicle development remains largely elusive.


A total of 20,874 mRNAs, 25,831 circRNAs, 4087 lncRNAs, and 794 miRNAs were annotated. Among them, we identified 58 DE lncRNAs, 325 DE circRNAs, 924 DE mRNAs, and 228 DE miRNAs during the development of medullary primary hair follicle development. GO and KEGG functional enrichment analyses revealed that the JAK-STAT, TGF-β, Hedgehog, PPAR, cGMP-PKG signaling pathway play crucial roles in regulating fibroblast and epithelial development during skin and hair follicle induction. Furthermore, the interactive network analysis additionally identified several crucial mRNA, circRNA, and lncRNA molecules associated with the process of primary hair follicle development. Ultimately, by investigating DEmir’s role in the ceRNA regulatory network mechanism, we identified 113 circRNA–miRNA pairs and 14 miRNA–mRNA pairs, including IGF2BP1-miR-23-x-novel-circ-01998-MSTRG.7111.3, DPT-miR-370-y-novel-circ-005802-MSTRG.14857.1 and TSPEAR-oar-miR-370-3p-novel-circ-005802- MSTRG.10527.1.


Our study offers novel insights into the distinct expression patterns of various transcription types during hair follicle morphogenesis, establishing a solid foundation for unraveling the molecular mechanisms that drive hair development and providing a scientific basis for selectively breeding desirable wool-related traits in this specific breed.

Peer Review reports


The wool fibers of Tibetan sheep, derived from primary hair follicles, primarily consist of myeloid wool, making carpet wool a valuable genetic resource known for its high quality. There is currently widespread interest in the study of hair follicle development in sheep and awareness of its importance for understanding the adult fleece and controlling wool production [1]. The initial stage of hair follicle development is a complex process, commencing with the proliferation of epidermal cells to form a placode beneath which an aggregation of dermal cells occurs, and these two cell populations interact and undergo orderly proliferation and differentiation. Subsequently, clumps gradually emerge within the dermis, ultimately giving rise to primary hair follicles, while certain newly formed hair follicles progress into secondary ones [2]. The rapidly proliferating matrix cells in the hair bulb produce hair shafts, and the dermal papilla consists of specialized fibroblasts located at the base of the hair follicle, which control the number of matrix cells and thus the size of the hair [3]. Primary hair follicle is characterized by the presence of large sebaceous glands and sweat glands, with both the diameter of hair follicles and hair bulbs being larger than those of secondary hair follicles, which develop earlier. However, all types of hair follicles possess complete sebaceous glands that contain a certain proportion of hollow or myelinated rough and elastic hairs [2, 4]. The morphology structure and development of hair follicles have important effects on the economic performance of wool. Studying the molecular mechanisms of sheep hair follicle morphogenesis can provide a refined and highly instructive research model for addressing key problems in modern biology [2], and will also contribute to understanding the genetic basis of wool related traits.

As the complexity of trait analysis increases, non-protein-coding sequences are increasingly prevalent in the genomes of multicellular organisms. Various classes of small and large non-coding RNAs (ncRNAs) have been demonstrated to regulate gene expression at almost every level, encompassing chromosome dynamics, chromatin architecture, transcriptional control, post-transcriptional processing, and translation, focusing primarily on animals [5]. Research over the last two decades revealed new classes of ncRNAs, including microRNAs (miRNAs), long ncRNAs (lncRNAs), and circular RNAs (circRNAs), all with different regulatory functions that are interwoven together in the larger RNA communication network to regulate the fundamental protein effectors of cellular function [6], including the regulation of hair follicle morphogenesis and development through DNA-binding transcriptional regulators, transcriptional activation, RNA editing, RNA interference and translational regulation [7,8,9,10]. Recent studies have shown that RNAs influence each other’s levels by competing for a limited pool of microRNAs, giving rise to a new theory of “competitive endogenous RNA” (ceRNA). mRNAs, lncRNAs, miRNAs, and circRNA activities form a large-scale regulatory network across the transcriptome, greatly expanding the functional genetic information in the genome and playing an important role [11]. Genome-wide identification of mRNAs, lncRNAs, circRNAs and miRNAs can take significant effect in processes of bovine intra-muscular adipogenesis [12]. circRNA-miRNA-mRNA ceRNA regulatory network affects muscle growth and development and individual growth by regulating the transformation of oxidized muscle fibers and glycolytic muscle fibers in Tibetan sheep [13].

Several key signals in WNT, BMP, EDAR, and FGF pathways, and a series of lncRNAs, were shown to be potentially important roles of in primary hair follicle induction and skin development [14]. More and more new evidence suggests the existence of a large network of regulatory ncRNAs, which constitute the sequence programming regulatory information of biological organisms and play an important role in the epigenetic development of their traits. Moreover, the variation of complex characteristic sequences primarily resides in the non-coding region, presumably the regulatory region, which is likely to alter the interaction between regulatory RNAs and their target, a prospect needs to be taken into consideration [15]. However, there are few reports about the transformation mechanism of hair follicle and skin development underlying regulation of the ceRNA network in carpet sheep fetal skin at placode stages.

To investigate the key genes, signaling pathways, and regulatory networks involved in medullary wool morphogenesis, our focus was directed towards analyzing the gene expression profiles of primary and secondary hair follicles. We utilized RNA-seq technology to investigate the regulatory role of competing endogenous RNAs (ceRNA) in the development of hair follicles in sheep skin and constructed a ceRNA regulatory network to identify key factors involved in hair follicle formation. This study will serve as a foundation for further research on hair follicles, provide insights into the interaction networks of ceRNA during the initiation of primary hair follicles, and provide a scientific basis for the breeding of wool-related traits.

Materials and methods

Animals and samples preparation

Carpet wool Tibetan sheep from Qinghai sheep breeding and promotion service center were used in this study. The study selected six healthy Tibetan ewes (3 years old; mean fiber diameter of coarse wool, 58.98 ± 2.48 μm; cashmere wool, 18.71 ± 2.73 μm) that did not exhibit estrus during the second estrus period following artificial assisted mating with the same ram, and these ewes were maintained under identical conditions. Based on previous morphological studies on the genesis and development of hair follicles [2], fetal skin samples were collected in pregnant ewes at E65 (TWP) and E85 (TWS). These samples were immediately placed in liquid nitrogen (− 80 ℃) for the transcriptome sequencing and Real-time quantitative polymerase chain reaction analysis (RT-qPCR).

All skin tissue (diameter 2 cm) was collected from the right medial area behind the shoulder blade to make a section and rinsed with 1 x phosphate buffered saline (PBS), then prepared paraffin-embedded sections and fixed in 4% paraformaldehyde at 4 °C for approximately one week before H&E staining.

RNA extraction, library construction

Total RNA was extracted using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. The quality and quantity of RNA were verified with 1% agarose gel, a NanoPhotometer® spectrophotometer (Implen, USA), Qubit® 2.0 Fluorometer (Life Technologies, USA), and Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). A total of 1.5 µ g RNA was used for sequencing library construction. The RNA-seq libraries were constructed using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) utilizing the manufacturer’s instruction. After purification and end-repair, the poly-A tail and sequencing adapters were ligated. Ligation products were size-selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq4000 by Gene Denovo Biotechnology Co (Guangzhou, China).

circRNA, lncRNA and miRNA library construction and sequencing

Fastp [16] (version 0.18.0) was used to filter adapters or low quality bases from raw reads. Short reads alignment tool Bowtie2 [17] (version 2.2.8) was used for mapping reads to ribosome RNA (rRNA) database. The rRNA removed reads of each sample were then mapped to reference genome (ncbi_GCA_002742125.1_Oar_v1.0) by HISAT2 [18].

Identification of circRNA: 20mers from both ends of the unmapped reads were extracted and aligned to the reference genome to find unique anchor positions within splice site. Anchor reads that aligned in the reversed orientation (head-totail) indicated circ RNA splicing and then were subjected to find_circ to identify circ RNAs [19].

Identification of lncRNA: All reconstructed transcripts were aligned to reference genomes to obtain new transcripts. We used the following parameters to identify reliable novel genes: the length of transcript was longer than 200 bp and the exon number was more than two. Three softwares CNCI [20] (version 2) and CPC [21] (version 0.9-r2) and FEELNC [22] (version v0.2) were used to choose long non-coding RNAs.

Identification of miRNA: All of the clean tags were aligned with small RNAs in GeneBank database (Release 209.0) to identify and remove rRNA, scRNA, snoRNA, snRNA and tRNA. All of the clean tags were then searched against miRBase database (Release 22) to identify known miRNAs. The known miRNAs were identified by comparison with miRNAs alignment from other species. All of the unannotated tags were aligned with reference genome. According to their genome positions and hairpin structures predicted by parameters of software mirdeep2, the novel miRNA candidates were identified.

Differentially expressed RNAs

To identify differentially expressed (DE) lncRNAs, mRNAs, circRNAs and miRNAs across three biological replicates samples, the edgeR package software between two different groups was used [23]. The fragments per kilobase per million mapped fragments (FPKM) value was calculated to quantify expression abundance of lncRNA and mRNA, the reads per million mapped reads (RPM) were used to compute expression abundance of the identified circRNAs, and the miRNA expression level was calculated and normalized to transcripts per million (TPM). Significant differential expression lncRNA and mRNA were identified with |log2Ratio| ≥ 1 and qval < 0.05 as screening conditions for significant DE genes. For circRNAs, fold change ≥ 2 and pval < 0.05 were used as screening parameter for significant DE circRNAs. miRNAs were identified with pval ≤ 0.05, and fold change ≥ 2 in a comparison as significant DE miRNAs.

RT-qPCR analysis

Several differentially gene expressed level involved in primary hair follicle induction were selected and confirmed by RT-qPCR. The housekeeping gene GAPDH and actin were used to normalize the expression levels of the lncRNAs, circRNAs and mRNAs, and U6 was used to normalize the expression levels of miRNAs. The primer sequences of the selected lncRNAs, mRNAs, circRNAs and miRNAs are listed in Supplementary Table 1. Quantitative real-time PCR was carried out in a 20 µL reaction mixture containing 10 µL 2×Top Green EX-Taq Mix, 2 µL cDNA, 7 µL ddH2O, and 0.5 µL forward and reverse primers. The following thermocycling program was used: 1 cycle of 94 ℃ for 30 s; 42 cycles of 94 ℃ for 5 s, 61 ℃ for 35 s; 1 cycle of 95 ℃ for 10 s, 65 ℃ for 60 s, 97 ℃ for 1 s followed by a final melting curve analysis stage. The relative expression levels were calculated using the 2−ΔΔ CT method [24]. All experiments were performed in three biological replicates.

GO and KEGG functional enrichment analysis

Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were analyzed using the Gene Ontology database (, and KEGG was the major public pathway-related database ( The calculating formula of P-value is: \(\text{P}=1-\sum _{i=0}^{m-1}\frac{\left(\frac{M}{i}\right)\left(\frac{N-M}{n-i}\right)}{\left(\frac{N}{n}\right)}\), here N is the number of all genes with GO/KEGG annotation; n is the number of source genes in N; M is the number of all genes that are annotated to the certain GO terms / specific pathways; m is the number of source genes in M.

Interaction Network Construction

For prediction of mRNAs interacting with circRNAs, lncRNAs and miRNAs, miRTarBase [25] (version 6.1) was used to predict mRNAs targeted by miRNAs sponge. The resulting correlation of circRNAs-miRNAs-lncRNAs-mRNAs can be visualized by Cytoscape [26] software to present a core and hub gene biological interact. The construction of ceRNA regulatory network complies with the following three aspects : (1) co-expressed negatively mRNA-miRNA or lncRNA / circRNA-miRNA expression correlation pairs was chosen using the Spearman Rank correlation coefficient (SCC) < -0.7. (2) co-expressed lncRNA–mRNA pairs or circRNA-mRNA pairs were selected by the PCC > 0.9. (3) Hypergeometric cumulative distribution function test was selected to test common miRNA sponges between two genes with p values <0.05. The calculating formula of P-value is: \(\text{p v}\text{a}\text{l}\text{u}\text{e}=1-\text{F}\left(x/U,M,N\right)=1-\sum_{i=0}^{x-1}\frac{\left(\frac Mi\right)\left(\frac{U-M}{N-i}\right)}{\left(\frac UN\right)}\).


Histological characteristics of medullary hair follicles during development

Histological changes in early hair follicle development can be observed in horizontal and longitudinal skin sections stained with the haematoxylin-eosin (H&E). The formation of primary hair follicles (PFs) is indicated by the emergence of hair placode and dermal condensate at E65, while the sheep fetal skin epidermis appears thin and homogeneous (Fig. 1A, C and E). The formation of secondary follicles (SFs) commenced at E85, with PFs exhibiting earlier development. Additionally, PFs displayed increased size and length, while SFs appeared smaller and grew in close proximity to the PFs epidermis. Furthermore, the dermis of sheep skin exhibited greater thickness compared to the preceding period (Fig. 1B, D and F).

Fig. 1
figure 1

Hematoxylin-eosin (HE) staining of sheep hair follicles at two developmental stages. Horizontal (magnification: ×400) (A, B) and longitudinal (×400) (C, D) slices of skin of primary (A, C) and secondary hair follicles (B, D). Horizontal (×100) and longitudinal (×100) slices of skin of primary hair follicle (E, F). PF: Primary hair follicle; SF: Secondary hair follicle; Epi, epidermis

Overview of mRNAs/circRNAs/lncRNAs sequencing

After filtering out the low-quality reads, a total of 520,980,026 valid reads with Q30 greater than 91.46%were obtained in all samples, of which more than 88.99% were mapped to the reference genome of Ovis aries (Supplementary Table 1). A total of 20,874 mRNAs (novel-genes 373), 4087 lncRNAs (3315 existing and 772 novel), and 25,831 circRNAs with an average length of 891 bp as novel circRNAs were identified during medullary hair follicle development (Supplementary Table 2). Meanwhile, the identified lncRNAs contained 89.23% long intergenic noncoding RNA (lincRNAs), and 75.87% of typical circRNAs were identified (Fig. 2A, B). The length distribution of lncRNAs was consistent with that of protein-coding gene. Both lncRNAs and mRNAs contain exons from 1 to 20, most lncRNAs had 1–2 exons. (Fig. 2C). The distributions of lncRNAs and mRNAs lengths were consistent (Fig. 2D). The expression heatmaps of circRNAs, mRNAs, lncRNAs showed inter-group and intra-group differences and clustering differentiation (Fig. 2E, F, G). In addition, the violin plot showed differences in circRNA, mRNA, lncRNA expression distribution (Fig. 2H, I, J). Two groups showed a total of significantly differentially expressed 325 DE circRNAs (including significantly 173 up-regulated and 152 down-regulated), 299 circRNAs were derived from source genes, and residual 26 circRNAs were derived from intergenic regions, having no source genes, 924 DE mRNAs (including significantly 452 up-regulated and 472 down-regulated), and 58 DE lncRNAs (including significantly 20 up-regulated and 38 down-regulated) (Fig. 2K, L, M).

Fig. 2
figure 2

Expression patterns of mRNAs, circRNAs and lncRNAs. A The type and proportion of circRNAs. B The type and proportion of lncRNAs. C Exon number distribution of lncRNAs and mRNAs. D Transcript lengths distribution of lncRNAs and mRNAs. E-G. Clustered expression heatmaps of all mRNAs. H-J. Violin plot of gene expression levels all mRNAs. K-M. The Venn diagrams of the shared and unique differential expressed circRNAs, mRNAs, lncRNAs

Differential expression analysis of mRNA/circRNA/lncRNA

To discover the biological functions of differentially expressed genes (DEGs) during early hair follicle development, functional enrichment analysis of stage-specific DEMs, host genes of DECs and cis-regulated genes of DELs were performed. The results showed that the GO terms and KEGG signaling pathways of DEMs were mainly involved in PPAR signaling pathway, ECM-receptor interaction, Retinol metabolism, cAMP signaling pathway, biological adhesion, and cell adhesion (Fig. 3A and B, Supplementary Tables 3 and 4). In addition, we found that some pathways were significantly enriched in up-regulated and down-regulated genes, respectively (Fig. 3G and H, Supplementary Tables 5 and 6), after further GSEA analysis of these pathways by Dex treatment, we found that the Circadian entrainment and Glycosphingolipid biosynthesis pathways were significantly activated, the PPAR signaling pathways, ECM receptor interaction, and Retinol metabolism were suppressed (Fig. 3I). These pathways were highly related to the early stage of hair follicle induction. DECs might be mainly enriched in Phosphatidylinositol signaling system, Hedgehog signaling pathway, ECM-receptor interaction, clathrin coat and growth factor binding (Fig. 3C and D, Supplementary Tables 7 and 8). The results of functional enrichment analysis of DELs indicated that DELs might be involved in cGMP-PKG signaling pathway, Prolactin signaling pathway and phospholipid biosynthetic process (Fig. 3E and F, Supplementary Tables 9 and 10).

Fig. 3
figure 3

GO functional annotation and KEGG enrichment analysis. A-B. The top 20 KEGG and GO enrichment pathways for DE mRNAs. C-D. The top 20 KEGG and GO enrichment pathways for DE circRNAs. E-F. The top 20 KEGG and GO enrichment pathways for DE lncRNAs. G-H. Pathways enriched by up-regulated or down-regulated DE mRNAs. I. GSEA enrichment analysis of selected results

Construction of the DE circRNA/lncRNA–mRNA interaction network

To further elucidate the key genes involved in primary hair follicle development, we subsequently constructed a circRNA/lncRNA-mRNA network interaction map based on the identified shared differential genes (Fig. 4, Supplementary Table 11). A total of 64 potential DECs and 13 key DELs that are likely to play crucial roles in the development of primary hair follicles. The findings revealed that the interaction network comprised a total of 1056 connections, with each lncRNA exhibiting associations with multiple mRNAs. A total of 12 circRNAs and three lncRNAs were identified as co-expressed with Dsc1, while 13 circRNAs and one lncRNAs were found to exhibit co-expression with ALX4, the network also includes circRNAs/lncRNAs that are co-expressed with Egr3, TGM3, KRT77, HOXC13, TSPEAR, MEGF9, SHH, SOX9, PTGER2. The aforementioned circRNAs/lncRNAs associated with genes may exert a significant influence on the regulation of hair follicle formation and growth.

Fig. 4
figure 4

circRNA/lncRNA-mRNA interaction network. The edges are the Pearson correlation coefficient (PCC) between genes. DE mRNAs and DE circRNAs / lncRNAs are represented by ellipse and diamonds, respectively

Identification of a miRNA landscape

By removing other non-coding RNAs such as rRNA, scRNA and tRNA, a total of 794 miRNAs were obtained (Fig. 5A). The initial base of these miRNAs was uridine and they were mostly expressed at the middle (0.01 < percentage < 1) levels (Fig. 5B). Violin plot and clustered expression heatmap of miRNA showed expression distribution differences and clustering differentiation (Fig. 5C and D). There were 156 upregulated and 72 downregulated DE miRNA that were identified (Fig. 5E and F). The biological functions of the target genes of these DE miRNAs were further analyzed by functional enrichment. The GO terms and KEGG signaling pathways were mainly involved in ECM-receptor interaction, PI3K-Akt signaling pathway, Growth hormone synthesis, secretion and action, VEGF signaling pathway, metabolic process, intracellular organelle (Fig. 5G and H, Supplementary Tables 12 and 13).

Fig. 5
figure 5

The landscape of miRNA. A The tag type of miRNAs. B First base preference of different tag lengths. C Violin plot of miRNAs expression levels. D Clustered expression heatmap of all miRNAs. E Scatter plots of all miRNAs. Red represents up-expressed miRNAs. Green represents down- expressed miRNAs. Blue is equally-expressed miRNAs. F Venn diagram showing the up-regulated and down-regulated DE miRNAs. G-H. The top 20 KEGG and GO enrichment pathways analysis of DE miRNAs

Construction of the DE circRNA/lncRNA–miRNA–mRNA ceRNA regulatory network

The regulatory network of circRNA/lncRNA-mRNA interaction integrates the prediction results of miRNA target genes, thereby facilitating the exploration of regulatory elements governing primary hair follicles. The participation of DEmiRs in the ceRNA regulatory network was analyzed, and subsequently integrated with the circRNA/lncRNA-mRNA regulatory relationship to construct a comprehensive ceRNA network (Fig. 6, Supplementary Table 14). The ceRNA regulatory network contained 113 circRNA–miRNA pairs and 14 miRNA–mRNA pairs and included 6 mRNAs, 44 circRNAs, 24 lncRNAs and 9 miRNAs. MSTRG.7111.3, novel-circ-019989 and miR-23-x have shared target gene IGF2BP1, Similarly, we observed the same findings in DPT-miR-370-y-novel-circ-005802-MSTRG.14857.1 and TSPEAR-oar-miR-370-3p-novel-circ-005802-MSTRG.10527.1. The network regulation also involves HSD11B2, CNKSR2, Rarb, and miR-195-x. The ceRNA network may offer valuable insights into the induction of primary hair follicle formation.

Fig. 6
figure 6

The ceRNA co-regulation network. The circle, triangle, diamonds, and V-type represented the DE mRNAs, DE lncRNAs, DE circRNAs, and DE miRNAs. The higher the connectivity, the redder the color

Differential expression of ceRNA-RT-qPCR validation

The selection of four mRNAs (DPT, IGF2BP1, Rarb, and CNKSR2), four lncRNAs (ENSOART00020016315, ENSOART00020033476, MSTRG.14334.1 and MSTRG.8125.2), and four circRNAs (novel_circ_000329, novel_circ_009066, novel_circ_002228, and novel_circ_025367) and 4 miRNAs (miR-370-y, miR-322-x, miR-497-x, and miR-10,174-x) were conducted validate the sequencing results using real-time quantitative PCR (Fig. 7). The consistency between their expression trends and sequencing results suggests the robustness of the gene expression profile data obtained in this study.

Fig. 7
figure 7

The expression level of four transcript types (mRNA, circRNA, lncRNA and miRNA) in the ceRNA network were validated by RT-qPCR. The left y-axis represents the RT-qPCR data, while the right side displays the RNA-seq data


The wool hair follicle serves as a unique and invaluable model for fundamental research in the life sciences, offering fertile ground to explore the underlying mechanisms of hair follicle morphogenesis, regeneration, heredity related to hair growth, and hair structure. The previous studies have demonstrated that during the induction stage of primary hair follicles, the regulatory mechanisms governing wool and hair follicle induction are largely conserved; however, distinct morphological changes occur with different characteristic genes, including certain regulatory elements such as lncRNA, which exhibit non-conserved functions across species [27, 28].

Although the role of lncRNA in hair follicle and skin development has been investigated, the precise mechanism underlying hair follicle induction remains elusive. Nie et al. identified 62 differentially expressed lncRNAs in primary hair follicles and predicted their potential targets, including Wnt16, Bmp3, and Itgb4, which are significantly enriched in key pathways such as WNT signaling, Hedgehog signaling, focal adhesion, and ECM receptor interaction. These findings align with our own study results and further underscore the crucial involvement of these pathways in hair follicle morphogenesis and skin development [14]. LNC_002919 and novel_circ_0026326 were found to act as ceRNA and participated in the regulation of the HF cycle as miR-320-3p sponges during skin development [29]. LncRNA MSTRG14109.1 and circRNA452 were competed with miRNA-2330 to regulated the cashmere fineness [30]. The findings suggest that ncRNAs may have a multifaceted and pivotal role in the regulation of hair genesis and growth. Furthermore, conducting a comprehensive analysis involving multiple regulatory elements could enhance our understanding of this intricate biological process. In this study, our conducted a comprehensive analysis of full transcriptome sequencing across different stages of hair follicle development, with the objective of identifying the genetic components involved in primary hair follicle occurrence and development. Currently, there is a lack of research on ceRNA network regulation pertaining to wool fiber characteristics in Tibetan sheep from Qinghai Province. Therefore, we first constructed the ceRNA regulatory network within the skin tissue of coarse wool fetal sheep.

A total of 58 DE lncRNAs, 325 circRNAs, 924 mRNAs, and 228 miRNAs were identified through high-throughput sequencing, playing a pivotal role in facilitating the development of primary hair follicles. The number of circRNAs annotated in primary hair follicles in this study was comparatively lower than the number annotated in secondary hair follicle studies conducted on cashmere goats (876 circRNAs) [31], and higher than the non-coding RNAs sequencing results were obtained between small waves and straight wool groups (114 circRNAs) [32]. However, this result is similar to the results of systematic analysis of non-coding RNAs in Angora rabbits during three hair follicle (HF) growth stages (247 circRNAs) [29]. In terms of the types and proportions of ncRNA sequencing, long intergenic noncoding and typical circRNAs accounted for the largest proportion, which was consistent with previous studies on sheep fetal and postnatal hair follicle development and adheres to established patterns [33, 34].

JAK-STAT, TGF-β, Hedgehog, PPAR, cGMP-PKG signaling pathway have been identified as crucial regulators in development of fibroblasts and epithelium during skin and hair follicle induction as well as hair follicle morphogenesis and skin development hair follicle induction stimulation [35,36,37,38]. In our data, these pathways mentioned exhibit a significant enrichment of differentially expressed genes in both mRNA and ncRNAs, forming a complex genetic regulatory network. This suggests that these signaling pathways may play a crucial role in regulating the growth of hair follicles. Additionally, we observed a significant enrichment of DEMs, DECs and DELs in the Focal adhesion and the ECM receptor interaction pathways, mirroring the enrichment findings in goats [39, 40]. The WNT signaling pathway is active in both the epithelial and mesenchymal components of hair follicle development, and WNT signaling is essential for promoting proliferation of the follicular epithelium and facilitating dermal condensation into dermal papillae, maintaining the properties necessary for hair follicle induction playing a crucial role in controlling hair follicle development and circulation [41,42,43]. Our analysis revealed a significant enrichment of DEGs in the WNT pathway, highlighting the crucial role played by the WNT pathway in primary hair follicle dermal condensation formation. Additionally, the circadian entrainment and glycosphingolipid biosynthesis signaling pathways play a pivotal role in the regulation of hair follicle regeneration [44, 45].

To further elucidate the key candidate genes underlying primary hair follicle induction, we conducted additional gene identification, integrated DEGs with DELs and DECs, and constructed a comprehensive circRNA/lncRNA–mRNA cross-regulatory network map. The presence of Dsc1 has been observed between cells at various stages of terminal differentiation in hair follicle (HF) keratinocytes [46], and knockout mice demonstrate epigenetic-follicular degeneration in older individuals [47]. The ALX4 gene plays a crucial role in the development of hair follicles, and individuals with homozygous c.793 C→T nonsense mutations in this gene may experience complete alopecia [48]. The gene HOXC13 may exert regulatory effects on hair follicle circulation in the context of LNC_000881 [49] and have been pinpointed as key determinants of wool fineness in sheep [50]. The protein KRT77, belonging to the type II keratin family, plays a crucial role in the formation of epidermis and coat formation [51]. Competitive adsorption of miR-184 liberates FGF10, which is implicated in hair follicle development, from chi-circRNA-0001141 [52]. Network analysis revealed the presence of genes such as TSPEAR [53], Egr3 [54], MEGF9 [55], SHH [56], PTGER2 [57], TGM3 [58], S100A7 [59], and SOX9 [60] that are intricately associated with hair follicle development and wool-related traits. The findings of these studies suggest that the screened wool-related traits candidate mRNAs and ncRNAs may unveil key factors involved in hair follicle growth, thereby serving as potential targets for subsequent investigations.

The regulation of targeted mRNAs expression by miRNAs governs a diverse array of cellular processes and functions, while also engaging in interactions with ncRNAs to form an intricate and interconnected ceRNA network, thereby playing a pivotal role in the regulation of target genes [61]. The expression and activity of the identified miR-195 [62], miR-370-3p [63], miR-432 [64], the miR-200 family [65], and miR-21 [66] can be regulated by target gene products. Changes in RNA abundance have the potential to impact ceRNA interactions [67], ultimately influencing hair follicle formation and development. Consequently, these findings suggest that candidate genes associated with various wool fiber traits are implicated in these molecular events. Our findings unveil several potential ceRNA regulatory networks involved in the process of hair follicle development and fiber growth. It has been reported that IGF2BP1 [68], HSD11B2 [69], TSPEAR [53], DPT [70], CNKSR2 [66], Rarb [71], and miR-195-x [62] play crucial roles in regulating hair follicle growth. Furthermore, the expression levels of IGF2BP1 [31], CNKSR2 [66] and Rarb [72] are intricately regulated through competitive binding interactions between lncRNAs/circRNAs and miRNAs during the process of hair follicle formation or signaling pathway activation. The reliability of gene expression in the ceRNA network was also confirmed by our RT-qPCR analysis. Apparently, the interaction between lncRNA/circRNA-mRNA and miRNA, as well as the targeting relationship between miRNA and mRNA in the ceRNA network, unequivocally underscores the pivotal role of miRNA in primary hair follicle and skin development of sheep, highlighting the intricate nature of the mechanism governing hair follicle initiation and development. Further investigations are required to elucidate the underlying mechanisms of action of these candidate mRNAs and non-coding RNAs in hair follicle induction, as well as to incorporate additional omics data (such as DNA methylome, ATAC-Seq, single-cell sequencing, and spatial transcriptomes) for the identification of molecular drivers involved in the regulation of hair follicle growth. In conclusion, a more comprehensive genomic data will enable us to gain a better understanding of the genetic basis underlying complex trait variation and pave the way for selective breeding strategies targeting wool traits in this particular breed.


In summary, we employed whole-transcriptome sequencing technology for the first time to characterize the expression profiles of mRNAs, lncRNAs, circRNAs, and miRNAs during primary hair follicle development in. coarse sheep fetal skin. The present study unveils crucial molecular features underlying the primary hair follicle formation and underscores the intricate regulatory networks between coding and non-coding genes during hair follicle development. Additionally, it unveils the underlying molecular mechanisms governing the interplay between coding and non-coding RNAs, thereby providing a valuable resource for comprehending the biology of hair follicle development and unraveling the genetic foundation of wool-related traits. Moreover, it facilitates subsequent coupling of molecular element that drive their wool-related traits morphological manifestation.

Availability of data and materials

The raw reads produced in this study were deposited in the NCBI SRA with the accession number SRA SUB13911514 under Bio-project PRJNA1033028. Additional data can be found in supplementary files.



Non-coding RNA


Circular RNA


Long non-coding RNAs


Micro RNA


Circular intronic RNA


Competing endogenous RNAs


long intergenic noncoding RNA


fragments per kilobase per million mapped fragments


reads per million mapped reads


transcripts per million


Gene Ontology


Kyoto Encyclopedia of Genes and Genomes


Pearson correlation coefficient


Spearman Rank correlation coefficient


differentially expressed genes


  1. Hardy MH, Lyne AG. The pre-natal development of wool follicles in Merino sheep. Aust J Biol Sci. 1956;9:423–41.

    Article  Google Scholar 

  2. Rogers GE. Biology of the wool follicle: an excursion into a unique tissue interaction system waiting to be re-discovered. Exp Dermatol. 2006;15:931–49.

    Article  PubMed  Google Scholar 

  3. Paus R, Cotsarelis G. The biology of hair follicles. N Engl J Med. 1999;341:491–7.

    Article  CAS  PubMed  Google Scholar 

  4. Schinckel PG. Follicle development in the Australian Merino. Nature. 1953;171:310–1.

    Article  CAS  PubMed  Google Scholar 

  5. Amaral PP, Mattick JS. Noncoding RNA in development. Mamm Genome. 2008;19:454–92.

    Article  CAS  PubMed  Google Scholar 

  6. Adams BD, Parsons C, Walker L, Zhang WC, Slack FJ. Targeting noncoding RNAs in disease. J Clin Invest. 2017;127:761–71.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Adelman K, Egan E. More uses for genomic junk. Nature. 2017;543:183–5.

    Article  CAS  PubMed  Google Scholar 

  8. Sun L, Goff LA, Trapnell C, Alexander R, Lo KA, Hacisuleyman E, et al. Long noncoding RNAs regulate adipogenesis. Proc Natl Acad Sci. 2013;110:3387–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Faghihi MA, Wahlestedt C. Regulatory roles of natural antisense transcripts. Nat Rev Mol Cell Biol. 2009;10:637–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Bai Y, Dai X, Harrison AP, Chen M. RNA regulatory networks in animals and plants: a long noncoding RNA perspective. Brief Funct Genomics. 2015;14:91–101.

    Article  CAS  PubMed  Google Scholar 

  11. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA. Language? Cell. 2011;146:353–8.

    Article  CAS  PubMed  Google Scholar 

  12. Huang C, Ge F, Ma X, Dai R, Dingkao R, Zhaxi Z, et al. Comprehensive analysis of mRNA, lncRNA, circRNA, and miRNA expression profiles and their ceRNA networks in the Longissimus Dorsi muscle of cattle-yak and Yak. Front Genet. 2021;12:772557.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Bao G, Zhao F, Wang J, Liu X, Hu J, Shi B, et al. Characterization of the circRNA–miRNA–mRNA network to reveal the potential functional ceRNAs Associated with dynamic changes in the Meat Quality of the Longissimus Thoracis muscle in Tibetan Sheep at different growth stages. Front Vet Sci. 2022;9:803758.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Nie Y, Li S, Zheng X, Chen W, Li X, Liu Z, et al. Transcriptome reveals long non-coding RNAs and mRNAs involved in primary wool follicle induction in carpet sheep fetal skin. Front Physiol. 2018;9:446.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Mattick JS. The genetic signatures of noncoding RNAs. Plos Genet. 2009;5:e1000459.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495:333–8.

    Article  CAS  PubMed  Google Scholar 

  20. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41:e166–166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35 suppl2:W345–9.

    Article  Google Scholar 

  22. Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, et al. FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017;45:e57–57.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  24. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2 [– Delta Delta C [T]] normalized to glyceraldehyde-3-phosphate dehydrogenase levels. qRT-PCR was method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  25. Hsu S-D, Lin F-M, Wu W-Y, Liang C, Huang W-C, Chan W-L, et al. miRTarBase: a database curates experimentally validated microRNA–target interactions. Nucleic Acids Res. 2011;39 suppl1:D163–9.

    Article  Google Scholar 

  26. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Nie Y, Li S, Zheng XT, Chen W, Li X, Liu Z, et al. Transcriptome reveals long non-coding RNAs and mRNAs involved in primary wool follicle induction in carpet sheep fetal skin. Front Physiol. 2018;9:1–16.

    Article  Google Scholar 

  28. Rinn JL, Chang HY. Long noncoding RNAs: molecular modalities to organismal functions. Annu Rev Biochem. 2020;89:283–308.

    Article  CAS  PubMed  Google Scholar 

  29. Zhao B, Chen Y, Hu S, Yang N, Wang M, Liu M, et al. Systematic analysis of non-coding RNAs involved in the angora rabbit (Oryctolagus cuniculus) hair follicle cycle by RNA sequencing. Front Genet. 2019;10:407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Hui T, Zheng Y, Yue C, Wang Y, Bai Z, Sun J, et al. Screening of cashmere fineness-related genes and their ceRNA network construction in cashmere goats. Sci Rep. 2021;11:21977.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Yin R, Wang Y, Wang Z, Zhu Y, Cong Y, Wang W, et al. Discovery and molecular analysis of conserved circRNAs from cashmere goat reveal their integrated regulatory network and potential roles in secondary hair follicle. Electron J Biotechnol. 2019;41:37–47.

    Article  CAS  Google Scholar 

  32. Lv X, Chen W, Sun W, Hussain Z, Chen L, Wang S, et al. Expression profile analysis to identify circular RNA expression signatures in hair follicle of Hu sheep lambskin. Genomics. 2020;112:4454–62.

    Article  CAS  PubMed  Google Scholar 

  33. 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:8501.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Zhao R, Liu N, Han F, Li H, Liu J, Li L, et al. Identification and characterization of circRNAs in the skin during wool follicle development in Aohan fine wool sheep. BMC Genomics. 2020;21:1–14.

    Article  CAS  Google Scholar 

  35. Harel S, Higgins CA, Cerise JE, Dai Z, Chen JC, Clynes R, et al. Pharmacologic inhibition of JAK-STAT signaling promotes hair growth. Sci Adv. 2015;1:e1500973.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Foitzik K, Paus R, Doetschman T, Dotto GP. The TGF-β2 isoform is both a required and sufficient inducer of murine hair follicle morphogenesis. Dev Biol. 1999;212:278–89.

    Article  CAS  PubMed  Google Scholar 

  37. Pummila M, Fliniaux I, Jaatinen R, James MJ, Laurikkala J, Schneider P, et al. Ectodysplasin has a dual role in ectodermal organogenesis: inhibition of Bmp activity and induction of Shh expression. Development. 2007;134:117–25.

  38. Bai T, Liang B, Zhao Y, Han J, Pu Y, Wang C, et al. Transcriptome analysis reveals candidate genes regulating the skin and hair diversity of Xinji fine-wool sheep and Tan sheep. Agriculture. 2021;12:15.

    Article  Google Scholar 

  39. Xu T, Guo X, Wang H, Hao F, Du X, Gao X, et al. Differential gene expression analysis between anagen and telogen of Capra hircus skin based on the de novo assembled transcriptome sequence. Gene. 2013;520:30–8.

    Article  CAS  PubMed  Google Scholar 

  40. Zhu B, Xu T, Yuan J, Guo X, Liu D. Transcriptome sequencing reveals differences between primary and secondary hair follicle-derived dermal papilla cells of the Cashmere goat (Capra hircus). PLoS ONE. 2013;8:e76282.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. DasGupta R, Fuchs E. Multiple roles for activated LEF/TCF transcription complexes during hair follicle development and differentiation. Development. 1999;126:4557–68.

    Article  CAS  PubMed  Google Scholar 

  42. Millar SE. Molecular mechanisms regulating hair follicle development. J Invest Dermatol. 2002;118:216–25.

    Article  CAS  PubMed  Google Scholar 

  43. Kishimoto J, Burgeson RE, Morgan BA. Wnt signaling maintains the hair-inducing activity of the dermal papilla. Genes Dev. 2000;14:1181–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Niu Y, Wang Y, Chen H, Liu X, Liu J. Overview of the circadian clock in the hair follicle cycle. Biomolecules. 2023;13:1068.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Bedja D, Yan W, Lad V, Iocco D, Sivakumar N, Bandaru VVR, et al. Inhibition of glycosphingolipid synthesis reverses skin inflammation and hair loss in ApoE–/– mice fed western diet. Sci Rep. 2018;8:11463.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Donetti E, Boschini E, Cerini A, Selleri S, Rumio C, Barajon I. Desmocollin 1 expression and desmosomal remodeling during terminal differentiation of human anagen hair follicle: an electron microscopic study. Exp Dermatol. 2004;13:289–97.

    Article  CAS  PubMed  Google Scholar 

  47. Chidgey M, Brakebusch C, Gustafsson E, Cruchley A, Hail C, Kirk S, et al. Mice lacking desmocollin 1 show epidermal fragility accompanied by barrier defects and abnormal differentiation. J Cell Biol. 2001;155:821–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Kayserili H, Uz E, Niessen C, Vargel I, Alanay Y, Tuncbilek G, et al. ALX4 dysfunction disrupts craniofacial and epidermal development. Hum Mol Genet. 2009;18:4357–66.

    Article  CAS  PubMed  Google Scholar 

  49. Wang S, Ge W, Luo Z, Guo Y, Jiao B, Qu L, et al. Integrated analysis of coding genes and non-coding RNAs during hair follicle cycle of cashmere goat (Capra hircus). BMC Genomics. 2017;18:1–13.

    Article  Google Scholar 

  50. Awgulewitsch A. Hox in hair growth and development. Naturwissenschaften. 2003;90:193–211.

    Article  CAS  PubMed  Google Scholar 

  51. Zhang M, Wu D, Ahmed Z, Liu X, Chen J, Ma J, et al. The genetic secrets of adaptation: decoding the significance of the 30-bp insertion in the KRT77 gene for Chinese cattle. Anim Biotechnol. 2023;34:3847–54.

    CAS  PubMed  Google Scholar 

  52. Gao Y, Song W, Hao F, Duo L, Zhe X, Gao C, et al. Effect of fibroblast growth factor 10 and an interacting non-coding RNA on Secondary Hair Follicle Dermal Papilla Cells in Cashmere goats’ follicle development assessed by whole-transcriptome sequencing technology. Animals. 2023;13:2234.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Peled A, Sarig O, Samuelov L, Bertolini M, Ziv L, Weissglas-Volkov D, et al. Mutations in TSPEAR, encoding a regulator of notch signaling, affect tooth and hair follicle morphogenesis. PLoS Genet. 2016;12:e1006369.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Wang Y, Zheng Y, Guo D, Zhang X, Guo S, Hui T, et al. m6A methylation analysis of differentially expressed genes in skin tissues of coarse and fine type liaoning cashmere goats. Front Genet. 2020;10:499712.

    Article  Google Scholar 

  55. Hlollis DE, Chapman RE. Apoptosis in wool follicles during mouse epidermal growth factor (mEGF)-induced catagen regression. J Invest Dermatol. 1987;88:455–8.

    Article  Google Scholar 

  56. Cui C-Y, Kunisada M, Childress V, Michel M, Schlessinger D. Shh is required for Tabby hair follicle development. Cell Cycle. 2011;10:3379–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Geng R-Q, Yuan C, Chen Y-L. Molecular cloning and expression analysis of prostaglandin E receptor 2 gene in cashmere goat (Capra hircus) skin during hair follicle development. Anim Biotechnol. 2014;25:98–107.

    Article  CAS  PubMed  Google Scholar 

  58. 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:1924.

    Article  CAS  PubMed  Google Scholar 

  59. Wu J, Li Y, Gong H, Wu D, Li C, Liu B, et al. Molecular Basis of Maintaining Circannual Rhythm in the Skin of Cashmere Goat. bioRxiv. 2020;2020.04.04.023044.

  60. Morita R, Sanzen N, Sasaki H, Hayashi T, Umeda M, Yoshimura M, et al. Tracing the origin of hair follicle stem cells. Nature. 2021;594:547–52.

    Article  CAS  PubMed  Google Scholar 

  61. Yamamura S, Imai-Sumida M, Tanaka Y, Dahiya R. Interaction and cross-talk between non-coding RNAs. Cell Mol life Sci. 2018;75:467–84.

    Article  CAS  PubMed  Google Scholar 

  62. Hui T, Zheng Y, Yue C, Wang Y, Bai Z, Sun J, et al. Screening of cashmere fineness-related genes and their ceRNA network construction in cashmere goats. Sci Rep. 2021;11:1–15.

    Article  Google Scholar 

  63. Zhao R, Li J, Liu N, Li H, Liu L, Yang F, et al. Transcriptomic analysis reveals the involvement of lncRNA–miRNA–mRNA networks in hair follicle induction in Aohan fine wool sheep skin. Front Genet. 2020;11:590.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Bai WL, Dang YL, Yin RH, Jiang WQ, Wang ZY, Zhu YB, et al. Differential expression of microRNAs and their regulatory networks in skin tissue of Liaoning cashmere goat during hair follicle cycles. Anim Biotechnol. 2016;27:104–12.

    Article  CAS  PubMed  Google Scholar 

  65. Hoefert JE, Bjerke GA, Wang D, Yi R. The microRNA-200 family coordinately regulates cell adhesion and proliferation in hair morphogenesis. J Cell Biol. 2018;217:2185–204.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Zhai B, Zhang L, Wang C, Zhao Z, Zhang M, Li X. Identification of microRNA-21 target genes associated with hair follicle development in sheep. PeerJ. 2019;7:e7167.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39:1033–7.

    Article  CAS  PubMed  Google Scholar 

  68. Xiao-Yan HE. Difference of MicroRNA expression in the ear and back skin of Young Alpaca (Lama Acos). Chin J Biochem Mol Biol. 2010;26:1016–22.

    Google Scholar 

  69. Boix J, Carceller E, Sevilla LM, Marcos-Garcés V, Pérez P. The mineralocorticoid receptor plays a transient role in mouse skin development. Exp Dermatol. 2016;25:69–71.

  70. Aoi N, Inoue K, Kato H, Suga H, Higashino T, Eto H, et al. Clinically applicable transplantation procedure of dermal papilla cells for hair follicle regeneration. J Tissue Eng Regen Med. 2012;6:85–95.

    Article  PubMed  Google Scholar 

  71. Everts HB, Sundberg JP, King LE Jr, Ong DE. Immunolocalization of enzymes, binding proteins, and receptors sufficient for retinoic acid synthesis and signaling during the hair cycle. J Invest Dermatol. 2007;127:1593–604.

    Article  CAS  PubMed  Google Scholar 

  72. Vettermann O, Siegenthaler G, Winter H, Schweizer J. Retinoic acid signaling cascade in differentiating murine epidermal keratinocytes: alterations in papilloma-and carcinoma‐derived cell lines. Mol Carcinog. 1997;20:58–67.

    Article  CAS  PubMed  Google Scholar 

  73. Percie du Sert N, Hurst V, Ahluwalia A, Alam S, Avey MT, Baker M, et al. The ARRIVE guidelines 2.0: updated guidelines for reporting animal research. J Cereb Blood Flow Metab. 2020;40:1769–77.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


We greatly appreciate our collaborators for their assistance in sample collection.


This work was supported by the Natural Science Foundation of Qinghai Province (2023-ZJ-748). National Breeding Joint Research Project.

Author information

Authors and Affiliations



KZ designed the experiments and revised the manuscript. DT analysed the data and wrote the paper. Others prepared the RNA samples and analysed the data. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Kai Zhao.

Ethics declarations

Ethics approval and consent to participate

This study followed the recommendations of the “Regulations for the Management of Affairs Concerning Experimental Animals” (Ministry of Science and Technology, China, revised in June 2004). The study was approved by the Animal Care and Use Committees of the Northwest Institute of Plateau Biology, Chinese Academy of Sciences. The animals are not harmed during sample collection. The study was conducted in accordance with the ARRIVE Guidelines [73].

Consent for publication

Not applicable.

Competing interests

The authors declare 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

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 The Creative Commons Public Domain Dedication waiver ( 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

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Tian, D., Pei, Q., Jiang, H. et al. Comprehensive analysis of the expression profiles of mRNA, lncRNA, circRNA, and miRNA in primary hair follicles of coarse sheep fetal skin. BMC Genomics 25, 574 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: