Expression profile analysis of microRNAs during hair follicle development in the sheep foetus

ABSTRACT

MicroRNAs (miRNAs) regulate the development and growth cycle of hair follicles (HFs). The molecular mechanism by which miRNAs determine the development of HFs in the sheep foetus remains elusive. In this study, the expression profiles of miRNAs at 11 development periods (45, 55, 65, 75, 85, 95, 105, 115, 125, 135 and 145 d) in sheep foetus skin were analysed by high-throughput sequencing and bioinformatics analysis. A total of 72 conserved miRNAs, 44 novel miRNAs and 32 known miRNAs were significantly differentially expressed. qRT-PCR results for 18 miRNAs were consistent with the sequencing data. 85 d of foetal development was the starting point for secondary hair follicle (SF) development according to tissue morphology and cluster analysis. In SF development, the prolactin signalling pathway and platelet activation played important roles, and 10 miRNAs were potential candidate miRNAs in SF initiation.

Graphical abstract

During 85d fetal development, 19 miRNAs were significantly up-regulated and 18 miRNAs were significantly down-regulated in initiation of secondary hair follicle growth

KEYWORDS:

• Hair follicle
• microRNAs
• development
• high-throughput sequencing

Hair follicles (HFs) are small but complex and are composed of 8 unique cell populations that are derived from the ectoderm and mesoderm. Thickening and elongation of epidermal keratino`cytes (hair placodes) occur in HFs [1], and several studies have suggested that HFs consist of a hair bulb, inner root sheath (IRS), hair shaft, outer root sheath, and connective tissue sheath [2,3]. HFoffine-wool sheep is comprised ofprimary hair follicles (PFs) and secondary hair follicles (SFs) [4–6], PFs occur early, whereas SFs occur late [7]. The growth and development of SFs directly influence the quantity and quality of wool produced [8,9]. SFs re-differentiation is the most important factor affecting HF density, which can be effectively increased thus improving the fineness of wool. The function and gene expression characteristics of HFs are distinct [10]. Initiation is the first and the most important stage of HF morphogenesis. Dermal cells send the initial signals to induce the formation of hair buds in the epidermis, followed by the release of TGF-β, EDA/EDAR/NF-κB, Noggin/Lef-1, shh,BMP-2/4/7,β-catenin and FGF in the hair buds to induce dermal fibrogenic cells to form the dermal papillae that regulate the activity of keratinocytes. The cells of dermal papillae release a second signal to stimulate epithelial cell proliferation, which initiates differentiation and results in a complete HF structure [11–13]. However, the precise timing of these processes and the molecular mechanism underlying HF initiation remain elusive.

Previous studies have primarily focused on gene function, whilethe roles of noncoding RNAs (such as long noncoding RNA, circular RNA, and microRNAs) at various stages have received increasing attentions.MicroRNAs (miRNAs) are a type of genetic material with approximately 22 bases in length and are in a class of endogenous, small, single-stranded noncoding RNA molecules. They play important roles in diverse biological processes by regulating posttranscriptional gene expression via pairing to the 3ʹ-untranslated regions (3ʹUTRs) of target mRNAs [14,15] and have a broad range of biological functions. Many studies have demonstrated that miRNAs play important roles in the regulation of HF development [16]. Yue et al.analysedthe roles of long noncoding RNA and coding RNA expression in sheep (Ovisaries) skin during the initiation of SF development by RNA sequencing (RNA-Seq) analysisandfound that the PPAR pathway likely plays a significant role during the initiation of SF development [17]. Liu et al.studied areas of the skin that produce cashmere during anagen using Solexa sequencing and showed that the miRNAswereco-expressed in goats and sheep andlocated in the same region of the respective chromosomes, they may play essential roles in skin and follicle development [18].

The Merino breed of sheep is prized for its wool, which hasthe finest and softest sheep wool.The Chinese Merino sheep (Junken Type) exhibit a high yield of quality wool with a fibre diameter (FD) of17.0–19.0μm. This study aimed to elucidate the regulatory mechanismof miRNA on Chinese Merino SF development in foetal period, which would benefit the cultivation of sheep with high-quality fine wool.

Materials and methods

Ethics statement

This study was conducted in strict compliance with the Guidelines for the Care and Use of Experimental Animals established by the Ministry of Science and Technology of the People’s Republic of China (approval number:2006–398) and theAnimal Care Guidelines and Regulationsfor the Administration of Laboratory Animals (2013 Revision) of the State Council of the People’s Republic of China, and it was approved by the Biological Studies Animal Care and Use Committee of Shihezi University, Xinjiang Academy of Agricultural and Reclamation Science. All efforts were made to minimize animal suffering.

Experimental animals and sample collection

A total of 50 2.5-year-old multiparous ewes ofthe super-fine-wool strainof Chinese Merino sheep (Junken Type)from the sheep farm of the Xinjiang Academy of Agricultural and Reclamation Science were used in this study.The sheep had an average wool FD of 17.5 ± 0.5μm according to the fibre testing standards of China and the International Wool Textile Organization (IWTO). All the ewes were sourced from a single flock and fed under the same conditions, and the oestrus time of eachwas determined using test rams. According to the pedigree information, the ewes were mated with the same strain of ram (FD = 16.6μm).The day of breeding was designated as d 0, and the gestational age of the foetuses was accurately recorded. A total of 33 foetuses were carefully, surgically removed from uteruses at gestational ages of 45d (period 1), 55d (period 2), 65d (period 3), 75d (period 4), 85d (period 5), 95d (period 6), 105d (period 7), 115d (period 8), 125d (period 9), 135d (period 10) and 145d (period 11); there were three replicates foreach time point. The foetuses were washed in phosphate-buffered saline and exsanguinated, and 1 cm² skin strips were rapidly removed from both mid-sides of each foetusand snap frozen in liquid nitrogen for RNA extraction. The remainder of each foetus was placed in 4% paraformaldehyde for morphological observation.

Paraffin section fabrication

The skin tissues of the foetuses were treated according to the conventional paraffin-embedded section preparation method. After dressing, the fixed skin was washed with PBS, dehydrated in concentrated ethanol, made transparent with xylene, and embedded in paraffin, and the embedded wax block was fixed on a slicer and continuously sectioned at a thickness of 7μm. Sections were removed with a glass slide, dried at 37°C, H.E. stained with a neutral resin and observed under a microscope; representative images were captured.

Total RNA extraction and small RNA library preparation

Total RNA was isolated from foetal skins using a TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Equal volumes of total RNA were detected by 1% formaldehyde denaturing agarose gel electrophoresis. The bands were clearly visible.Total RNA purity was quantified using an ND 1000 Nanodrop Spectrophotometer (Thermo Scientific, Wilmington, USA), and RNA integrity was verified using an Agilent 2100 Bioanalyzerwithan RNA 6000 Nano Chip (Agilent Technologies, Santa Clara, USA). The OD₂₆₀/OD₂₈₀ ratios were greater than 1.8, and the RIN values were greater than 8. The initial total amount of RNA was accurately quantified with a Qubit RNA Assay Kit (Invitrogen, Carlsbad, USA).Ligated 3ʹ and 5ʹ adapters from a TruSeq Small RNA Library Prep Kit (Illumina, San Diego, CA, USA) were ligated at both ends of the small RNAs in two separated reactions using T4 RNA ligase 2 truncated with ATP (NEB, Beijing, China). The reverse transcription (RT) reaction was performed using SuperScript II Reverse Transcriptase, and the cDNA was amplified with PCR master mix according to the manufacturer’s instructions. The libraries were diluted before being runon a High-Sensitivity DNA Chip (Agilent Technologies, Santa Clara, USA), purified by 6% PAGE and quantified with a Qubit DNA HS Assay Kit (Invitrogen, Carlsbad, USA).Library quality was assessed using an Agilent 2100 Bioanalyzer High-Sensitivity DNA Chip and quantified using KAPA SYBR FAST Master Mix Universal 2X qPCR Master Mix anda DNA Quantification Standards and Primer Premix Kit (KAPA Biosystem, Massachusetts, USA).The libraries were adjusted to a concentration of 2nM for sequencing witha HiSeq2500 System (Illumina, San Diego,CA,USA) following the manufacturer’s protocol.

Small RNA high-throughput sequencing and processing of reads

Total RNA was isolated from the skins of 33 fine-wool sheep foetuses, and equal quantities of three samples from the same developmental period were pooled and prepared according to the Small RNA Sample Prep Kit protocol. Illumina small RNA deep sequencing was performed at the Beijing CapitalBioCo.Ltd. Total RNA at approximately 2nM was providedfor each sample set, and individual sequence reads with base quality scores were produced by Illumina sequencing. Adapter sequences were removed from both ends of the reads, and all identical sequences were counted and eliminated from the initial data set. The resulting set of unique sequences, with their associated read counts, was referred to using the sequence tags. The unique reads were mapped onto the O. ariesgenome from Ensembl (http://ensembl.org) using the programBowtie [19].

Differential expression analysis and hierarchical clustering of microRNAs

Similar to the credibility interval approaches for analysing SAGE data [20], we employed edger to identify miRNAs showing statistically significant differences in relative abundance (as reflected by the total count of individual sequence reads) among 11 small RNA libraries [21]. Finally, miRNAs with a P-value<0.05 were considered to be significantly different among the 11 libraries. Differential miRNA TPM (transcripts per million) values were used to representthe level of expression differentiation among HFs at different foetal ages.

Hierarchical clustering analysis was performed with MEV software v4.6, and differently coloured regions were used to represent different clustering groups. Development and expression patterns within the same group were similar, possibly indicating similar functions or participation in the same biological processes.

Prediction of target genes for differential expression of microRNAs and functional analysis

A search for miRNA target genes was then performed using apreviously described approach. All newly identified miRNA sequences were used to query the 3ʹUTRs of O. aries for potential target sequences using RNAhybrid software with default parameters and permitting anenergy threshold of −20 kcal/mol. To determine the function of differentially expressed miRNAs in HF development, we subsequently analysed 3,223 genes through the DAVID website (http://david.abcc.ncifcrf.gov/). Furthermore, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation were conducted using the DAVID website.The significance threshold was set at a P-value<0.05.

Validation of small RNA sequencing by real-time PCR

The expression patterns of miRNA variedamong the eleven small RNA libraries constructed from foetal skin samples taken duringHF development. Eighteen miRNAs were randomly selected to confirm differential expression by quantitative real-time RT-PCR (qRT-PCR). Qualitative assessments of the purified RNAs were determined utilizing an ImplenNanophotometer UV Spectrophotometer (Germany) and gel electrophoresis, and total RNA was then reverse-transcribed using anmiScript II RT Kit (Qiagen,Germany). First-strand cDNA was synthesized from total RNA (2 μg of each sample)according to the manufacturer’s instructions, and the subsequent cDNA was used as a template for qRT-PCR, which was performed with a QuantiFast SYBR Green PCR Kit (Qiagen, Germany) on a LightCycler 480II Detection System (Roche, Germany). The PCR amplification was conducted in the thermal cycler using 20 μL of reaction mixture composed of 10 μL of 2x QuantiTect SYBR Green PCR Master Mix, 2 μL of 10x miScript Universal Primer, 2 μL of 10x miScript Primer Assay, and 2 μg of cDNA template. The PCR conditions were as follows: initial denaturation at 95°C for 5 min followed by 40 cycles at 98°C for 10 s, annealing at 57°C for 5 s (different miRNAs were slightly different), and extension at 72°C for 20 s. A melting curve was generated to confirm the identity of each PCR product. The entire experiment was performed three times for each developmental period, andeach sample was analysed three times. Expression was normalized to that of 5S, an internal housekeeping control. The primers used in theqRT-PCR are shown in Table 1. The miRNA downstreamprimer was the same universal primer as in the QuantiFast SYBR Green PCR kit(Qiagen,Germany); therefore, the table only shows the upstream primers.

Table 1. The primers of real-time PCR.

miRNA name	Primer sequences(5ʹ-3ʹ)	Annealing temperature (℃)
miR-184	TGGACGGAGAACTGATAAGGGT	55
oar-miR-381-3p	ATATACAAGGGCAAGCTCTCT	55
oar-miR-26a	TTCAAGTAATCCAGGATAGGCT	56
novel-14	TTTTGCAATATGTTCCTGAATAT	57
oar-miR-377-3p	ATCACACAAAGGCAACTTTCGT	57
oar-miR-431	TGTCTTGCAGGCCGTCATGCAGG	57
oar-miR-376a-5p	TAGATTCTCCTTCTATGAGTAC	57
oar-miR-133	TTGGTCCCCTTCAACCAGCTGT	57
miR-135b	TATGGCTTTTTATTCCTATGTGA	57
miR-504	AGACCCTGGTCTGCACTCTGTCT	57
miR-339b	TCCCTGTCCTCCAGGAGCTC	57
miR-206	TGGAATGTAAGGAAGTGTGTGG	57
miR-767	TGCACCATGGTTGTCTGAGCATG	56
oar-miR-323b	CCCAATACACGGTCGATCTCT	57
oar-miR-655-3p	ATAATACATGGTTAACCTCTCT	57
miR-208a	ATAAGACGAACAAAAGGTTTGT	57
miR-338*	AACAATATCCTGGTGCTGAGT	57
Novel-118	CACCCGTCCCGTGCGTCCCCGGA	57
5S-F	CTCGTCTGATCTCGGAAGCTAA	57
5S-R	CTACAGCACCCGGTATTC	57

Statistical analyses

For gene expression analysis, the number of unique-match reads was calculated and normalized to TPM values, and 18 differentially expressed miRNAs were selected according to the sequencing results. Gene expression was calculated and normalized to the 5Sexpression level, and the relative expression level of the target miRNA was calculated using the 2^−ΔΔCt method. The statistical software package SPSS17.0 was used to calculate the Ct values and standard errors among replicate samples. The expression levels of miRNAs in different developmental periods were compared using ANOVA.

Results

Morphological changesin hair follicles

The development of HFs was investigated by observing tissue sections. During the process of foetalHF development, epithelial cells and dermal cells initially interacted to form the primitive HF body, followed by the gradual growth and development of the complete HF structure.At 75d, the hair germ cells (a1) and dermal fibroblasts (a2) appeared in the foetal skin (Figure 1(a)). At 85d, SFs began to appear, and SFsfurther developed from the enlarged portion of the epidermal layer of the PFs (Figure 1(b)).

Figure 1. Variation in process of secondary hair follicle morphogenesis in Chinese Merino sheep during the foetal period.

A: 75 days of foetal development (100×); B: 85 days of foetal development (100×) a1: Hair germ cells, a2: Dermal fibroblasts, PF: Primary hair follicle, SF: Secondary hair follicle

Overview of the Solexa sequencing data

All sequencing data for each small RNA library are presented in Table 2. We mapped these clean reads using the sheep (O. aries) genome from Ensembl(http://ensembl.org) to further elucidate the possible mechanisms underlying the diversity of small RNAs involved in HF development. Mapping the total reads to parts of the genomeresulted in 77.4% of the clean reads in the eleven libraries (Table 3), and the distribution of small RNA length differedamong the eleven libraries. In all libraries, most of the clean reads (78.03%) were 22 nucleotides long (Figure 2).

Table 2. Small RNA sequence statistics for 11 libraries.

library		Total Reads	Low Quality	Adaptor3 Null	Insert Null	5ʹ Adaptor Contaminants	Size < 18 nt	PolyA	High Quality (size ≥ 18 nt)
1	#	10,630,196	16,202	41,929	4109	5890	308,221	68	10,253,777
	%	100	0.15	0.39	0.04	0.06	2.9	0	96.46
2	#	11,658,550	23,638	44,084	6044	4987	214,894	240	11,364,663
	%	100	0.2	0.38	0.05	0.04	1.84	0	97.48
3	#	11,406,751	24,444	42,193	4510	6480	288,645	139	11,040,340
	%	100	0.21	0.37	0.04	0.06	2.53	0	96.79
4	#	13,750,560	33,452	46,841	14,930	29,750	1,448,263	136	12,177,188
	%	100	0.24	0.34	0.11	0.22	10.53	0	88.56
5	#	15,653,293	49,273	5,037,571	6307	15,249	652,759	511	9,891,623
	%	100	0.31	32.18	0.04	0.1	4.17	0	63.19
6	#	16,046,032	51,398	5,829,044	5373	12,895	674,787	597	9,471,938
	%	100	0.32	36.33	0.03	0.08	4.21	0	59.03
7	#	10,138,981	23,446	44,540	3693	15,957	811,932	154	9,239,259
	%	100	0.23	0.44	0.04	0.16	8.01	0	91.13
8	#	10,146,459	20,196	49,614	5058	11,342	592,331	216	9,467,702
	%	100	0.2	0.49	0.05	0.11	5.84	0	93.31
9	#	23,285,212	46,993	83,358	16,997	23,231	1,426,534	257	21,687,842
	%	100	0.2	0.36	0.07	0.1	6.13	0	93.14
10	#	10,239,538	22,414	43,840	11,103	11,639	808,000	192	9,342,350
	%	100	0.22	0.43	0.11	0.11	7.89	0	91.24
11	#	10,552,874	18,803	40,860	7892	22,723	930,287	163	9,532,146
	%	100	0.18	0.39	0.07	0.22	8.82	0	90.33

#: the total number of reads; %: the percentage of the total number of reads.

Table 3. Expression Profile of Sequenced Reads in 11 libraries.

Class	High-quality Reads	Total Distinct Reads	Redundancy(%)(a)	Perfect Matched(b) Genome	%	Total of Perfect Matched	%
1	10,253,777	256,728	97.5	117,867	45.91	8,467,029	82.57
2	11,364,663	416,867	96.33	210,149	50.41	9,091,671	80
3	11,040,340	302,546	97.26	148,877	49.21	9,038,906	81.87
4	12,177,188	379,019	96.89	208,134	54.91	9,857,585	80.95
5	9,891,623	935,019	90.55	113,090	12.09	7,442,537	75.24
6	9,471,938	1,119,110	88.18	155,328	13.88	6,876,642	72.6
7	9,239,259	451,476	95.11	240,963	53.37	7,023,133	76.01
8	9,467,702	526,411	94.44	289,553	55.01	7,143,730	75.45
9	21,687,842	739,072	96.59	416,120	56.3	17,061,894	78.67
10	9,342,350	564,673	93.96	311,025	55.08	6,950,896	74.4
11	9,532,146	728,701	92.36	193,987	26.62	6,611,221	69.36
Total	123,000,000	4,841,912	96.08	1,720,096	35.53	95,565,244	77.4

Redundancy = 100-(Total Distinct Clean Reads/Total Clean Reads *100)

Figure 2. The distribution of the read lengths from the Solexa sequencing data.

Standard bioinformatics analysis was used to annotate the clean reads into different categories, including miRNA, known miRNA(69.1%), rRNA(2.5%), tRNA(4.0%), snoRNA(0.2%), and other RNA(0.4%). To identify the known miRNAs, unique sequences were aligned with all animal precursor and mature miRNA sequences from miRBase 18.0 (http://mirbase.org/), and the conservation of known miRNAs wasthen compared with those ofthe other species [22].

Identification of differentially expressed miRNAsamong different hair follicle development phases

The biogenesis of miRNAs is highly regulated during various time periods. Thus, comparing differentially expressed miRNAsamong libraries provides insight into the molecular events related to the developmental periods of HFs. We first normalized the readsand then used aP-value<0.05 and an absolute log2 ratio value ≥1.0 as thresholds to determine the statistical significance of differences in miRNA expression.The miRNAs that satisfied the criteria were denoted as differentially expressed miRNAs.Wehave listed the top ten most abundantly expressedknown miRNAs (Table 4), conserved miRNAs (Table 5) and novel miRNAs (Table 6).Among these miRNAs, oar-miR-10b had the highest expression ofthe known miRNAs, and miR-27 and Novel-14 had the highest expression ofthe conserved and novel miRNAs, respectively.

Table 4. The 10 most abundantly expressed known miRNAs in the skin and hair follicles.

microRNA	1	2	3	4	5	6	7	8	9	10	11
oar-miR-10b	350,828	352,254	317,900	263,991	236,946	188,878	182,164	197,927	175,706	172,212	185,980
oar-miR-26a	122,211	61,714	78,511	55,802	115,552	108,375	98,760	92,197	125,956	122,633	93,835
oar-miR-127	70,887	79,349	79,947	85,904	65,920	61,080	61,314	33,134	37,530	18,566	17,779
oar-miR-10a	58,146	38,240	45,463	23,207	24,584	20,907	23,683	27,516	47,218	46,892	32,145
oar-miR-148a	27,427	32,563	51,469	47,720	44,812	47,541	37,536	32,086	25,922	23,328	26,549
oar-miR-181a	22,695	22,997	29,128	38,530	55,547	62,245	46,627	57,781	50,454	43,995	43,142
oar-miR-143	22,121	19,056	23,506	41,635	44,321	38,482	45,212	44,688	59,842	53,202	97,345
oar-miR-411a-5p	16,703	23,551	27,257	14,102	21,160	20,913	20,094	10,809	11,472	5136	5299
oar-miR-199a-3p	16,283	14,000	17,334	25,220	20,462	22,672	14,505	12,856	13,777	11,933	13,517
oar-miR-30a-5p	9223	7690	7970	9090	6566	6431	8443	7872	7500	7514	9629

Table 5. The 10 most abundantly expressed conversed miRNAs in the skin and hair follicles.

microRNA	1	2	3	4	5	6	7	8	9	10	11
miR-27	5181	6736	10,024	12,900	21,166	27,838	25,536	26,932	32,034	33,240	40,792
miR-205	2009	1690	2250	4342	3624	3962	4166	4827	5154	5875	7349
miR-92	5118	5472	5735	9472	7790	8259	6080	6164	4466	5260	6704
miR-186	3456	4667	4332	3985	3752	4624	4176	5069	3504	4480	5230
miR-8	730	758	1072	1192	1917	2138	2167	2282	2566	3725	5035
miR-125b	1408	2978	2650	2500	3119	4051	2663	2771	3109	3276	3962
miR-151b	1337	1366	2315	1795	3837	4250	2530	2200	2894	2058	3060
miR-16	1727	1551	2892	2379	3846	3179	3943	3894	2773	2795	3056
miR-24	185	174	217	524	502	527	694	923	1393	1609	1822
miR-192	994	659	946	1512	1509	1339	1201	1393	1214	1169	1474

Table 6. The 10 most abundantly expressed novel miRNA in the skin and hair follicles.

microRNA	1	2	3	4	5	6	7	8	9	10	11
Novel-14	1124	2081	1745	974	1252	1117	761	749	453	308	268
Novel-125	122	131	104	137	102	84	62	67	51	43	42
Novel-117	36	49	48	153	63	100	171	167	167	212	132
Novel-124	28	53	50	131	59	121	132	100	103	149	59
Novel-31	27	41	40	30	35	33	34	31	26	18	21
Novel-38	25	27	31	28	29	31	22	31	19	18	17
Novel-110	24	14	6	7	11	10	6	7	10	6	3
Novel-9	21	24	14	40	16	31	66	72	63	96	47
Novel-25a	19	16	16	32	17	34	51	45	46	63	37
Novel-114	15	20	12	5	7	8	10	5	7	4	7

For the critical SF developmentperiod from 75–85 d, wehave listed the top ten most differentially expressed known miRNAs (Table 7), conserved miRNAs (Table 8) and novel miRNAs (Table 9). Six hundred seventy-nine (87.7%) miRNAswere co-expressed from 75–85 d; 61 (7.9%) were expressed at 75d; and 34 (85%) were expressed at 85d. For the differentially expressed miRNAs, 138 (95.8%) were co-expressed from 75–85 d; 4 (2.8%) were onlyexpressed at 75d; and 2 were expressed at 85d (1.4%). Four miRNAs were exclusively differentially expressed in the 75-d-library, namely, miR-146a, miR-196b*, miR-9*-1, and Novel-118*, and 2 miRNAs were exclusively differentially expressed in the 85-d-library, namely, miR-122 and miR-2285ak. A total of 37 miRNAs in sheep skin were differentially expressed between 75 and 85d of SF morphogenesis, among whichthe expression of 19 miRNAs was upregulated and that of 18 miRNAs was downregulated (Figure 3). Novel-126 showed the most significant downregulation in expression during SF development, whereas miR-196b showed the most significant upregulation.

Table 7. The 10 most abundantly expressed known miRNA in 75d and 85d of sheep hair follicles development.

microRNA name	75d	85d	Mature sequence	P value
oar-miR-10b	2,590,219	1,750,254	ACCCTGTAGAACCGAATTTGTG	0.97
oar-miR-127	842,871	486,929	ATCGGATCCGTCTGAGCTTGGCT	0.76
oar-miR-26a	547,518	853,553	TTCAAGTAATCCAGGATAGGCT	0.06
oar-miR-148a	468,218	331,016	TCAGTGCACTACAGAACTTTGT	0.89
oar-miR-143	408,513	327,389	TGAGATGAAGCACTGTAGCTC	0.68
oar-miR-181a	378,045	410,313	AACATTCAACGCTGTCGGTGAGT	0.28
oar-miR-199a-3p	247,457	151,150	ACAGTAGTCTGCACATTGGTT	0.85
oar-miR-10a	227,705	181,592	TACCCTGTAGATCCGAATTTG	0.69
oar-miR-411a-5p	138,368	156,300	ATAGTAGACCGTATAGCGTACG	0.24
oar-let-7i	132,557	54,407	TGAGGTAGTAGTTTGTGCTGTT	0.29

Table 8. The 10 most abundantly expressed conversed miRNA in 75d and 85d of sheep hair follicles development.

microRNA name	75d	85d	Mature sequence	P value
miR-27	126,568	156,350	TTCACAGTGGCTAAGTTCTGC	0.25
miR-92	92,934	57,541	TATTGCACTTGTCCCGGCCTGT	0.71
miR-205	42,600	26,771	TCCTTCATTCCACCGGAGTCTG	0.73
miR-186	39,103	27,718	CAAAGAATTCTCCTTTTGGGCT	0.94
miR-199	27,483	19,372	CCCAGTGTTTAGACTATCTGTTC	0.93
miR-125b	24,534	23,037	TCCCTGAGACCCTTTAACCTGT	0.59
miR-16	23,344	28,412	TAGCAGCACGTAAATATTGGAG	0.27
miR-100	19,218	14,220	AACCCGTAGATCCGAACTTGTG	0.98
miR-378a	18,638	9059	ACTGGACTTGGAGTCAGAAGG	0.37
miR-151b	17,613	28,346	TCGAGGAGCTCACAGTCTAGT	0.09

Table 9. The 10 most abundantly expressed novel miRNA in 75d and 85d of sheep hair follicles development.

microRNA name	75d	85d	Mature sequence	P value
Novel-14	9559	9251	TTTTGCAATATGTTCCTGAATAT	0.40
Novel-117	1505	467	TCTGGACGCGAGCCGGGCCC	0.10
Novel-125	1341	755	AGGGTTGGGCGGAGGCTTTCCT	0.73
Novel-124	1287	434	TTTCCCGGCCAATGCACCA	0.15
Novel-9	393	119	CGCGGGAGCCTCGGTCGGCCTCG	0.10
Novel-25a	315	128	GGCGCTGGAGCGTCGGGCCC	0.31
Novel-33	310	101	CCTCGACACAAGGGTTTG	0.14
Novel-31	297	255	ATGACTGTAATTCGCGTTAGCAGC	0.57
Novel-38	273	216	TCTCTGGGCCTGTGTCTTAGGCT	0.70
Novel-111	242	139	ACCCCACTCCTGGTACCA	0.78

Figure 3. Differentially expressed miRNAs in sheep skin between 75 and 85 d of secondary HF morphogenesis.

Hierarchical cluster analysis

To better understand the expression patterns of miRNAs during HF development, a cluster analysis was performed using MEV software v4.6 (Figure 4), and thedifferentially expressed miRNAs in 11 developmental periods were clustered into 2 categories and 4 subcategories.The development of PFs and SFs occurred from 45-95d.During this period, Primary hair follicle primordia was formed at first, then secondary hair follicles began to develop with primary hair follicles developing into the dermis at the embryonic age of 85 days Based on the foetuses at 95d, HFs continued to extend into the deep dermis withincreased age, gradually becoming longer and thicker. Theother categorywas 105-145d, during which mature sweat glands formed, most of the PFs and some of the SFs matured, and HFs formed keratinized hair shafts that penetrated the body surface. The four subcategories were as follows:the firstincluded 45, 55, and 65d of development, with 65dbeing the critical period for PF development;the second included 75, 85, and 95d of development, with 85d being the critical period for SF development;the third included 105, 115, and 125d of development; and the fourth included 135 and 145d of development.Overall, the cluster analysis results were consistentwith the critical stages for HF growth and development, which showedthat ourestablishedexperimental periods are meaningful.

Figure 4. Cluster analysis of the miRNA expression levels of during the 11 phases.

Predictive analysis of target genes of differentially expressed miRNAs

To understand the functions of differentially expressed miRNAs in hair follicle development, we predicted the target genes of 148 differentially expressed miRNAs using RNAhybrid predictive software, and the results showed that 932 genes were targeted by the differentially expressed miRNAs in the 11 stages ofHF development (Supplementary Table S1).The target genes of 37differentially expressed miRNAs were predicted using RNAhybrid to annotate their biological functions during SFdevelopment at 75 and 85d, and253 target genes were found to functionin SF development (Supplementary Table S2).

Functional analysis of target genes

To better understand the biological process of HF and skin morphogenesis, 582 target genes of differentially expressed miRNAs were categorized into three GO categories: cellular component, biological process and molecular function. GO analysis was conducted on differentially expressed genes in the 11 developmental periods of sheep skin HF morphogenesis, and the results showed that 33 genes in the cellular component, biological process and molecular function categories had significant differences in expression during the 11 periods (Table 10). Enrichment analysis found 19 genes that were differentially expressed in the cellular component, biological process and molecular function categories at 75 and 85 d of SF morphogenesis (Table 11). There were 5 common GO terms in the 11 SF development periods, including positive regulation of receptor internalization, nucleoplasm, nucleotide binding, poly(A) RNA binding, and ATP binding.

Table 10. Gene ontology analysis of differentially expressed genes in 11 periods development of sheep skin HF morphogenesis.

miRNA category	miRNA	GO_id	GO_term	Count	Go category	P-Value
Known miRNA	oar-miR-377-3p;oar-miR-3958-5p;oar-miR-376a-3p;oar-miR-539-3p.	GO:0060291	long-term synaptic potentiation	4	biological process	0.0029
Known miRNA	oar-miR-381-5p;oar-miR-544-3p;oar-miR-376c-5p	GO:0010875	positive regulation of cholesterol efflux	3	biological process	0.0089
Known miRNA	oar-miR-379-3p;oar-miR-411a-3p;oar-miR-381-5p;oar-miR-381-5p;oar-miR-377-3p; oar-miR-409-5p	GO:0008283	cell proliferation	6	biological process	0.0099
Known miRNA	oar-miR-133;oar-miR-133;oar-miR-381-5p;oar-miR-376c-5p	GO:0050728	negative regulation of inflammatory response	4	biological process	0.0202
Known miRNA	oar-miR-410-3p;oar-miR-23b	GO:0030718	germ-line stem cell population maintenance	2	biological process	0.0331
Known miRNA	oar-miR-376c-3p;oar-miR-411b-3p;oar-miR-377-3p;oar-miR-3958-5p	GO:0016192	vesicle-mediated transport	4	biological process	0.0418
Known miRNA	oar-miR-3956-5p;oar-miR-376c-5p	GO:0043031	negative regulation of macrophage activation	2	biological process	0.0439
Conserved miRNA	miR-504;miR-192;miR-215;miR-491*;miR-32*;miR-7c	GO:0042149	cellular response to glucose starvation	5	biological process	0.0015
Conserved miRNA	miR-9*-1;miR-345a;miR-133*;miR-196b*;miR-204a*;miR-7134*;miR-192;miR-215;miR-491*;miR-383	GO:0048147	negative regulation of fibroblast proliferation	5	biological process	0.0019
Conserved miRNA	miR-331*;miR-124;miR-450*;miR-92*; miR-32*;	GO:2,000,352	negative regulation of endothelial cell apoptotic process	4	biological process	0.0163
Conserved miRNA	miR-345a;miR-345b;miR-7857; miR-208a;miR-762a;miR-762b;miR-423*;miR-4454;	GO:0007268	chemical synaptic transmission	5	biological process	0.0230
Conserved miRNA	miR-345a;miR-345b;miR-192;miR-215;miR-491*;miR-383;	GO:0001836	release of cytochrome c from mitochondria	3	biological process	0.0409
Novel miRNA	Novel-118*;Novel-18a;Novel-18b; Novel-20	GO:0002092	positive regulation of receptor internalization	4	biological process	0.0018
Novel miRNA	Novel-118*;Novel-112;Novel-112	GO:2,000,679	positive regulation of transcription regulatory region DNA binding	3	biological process	0.0114
Novel miRNA	Novel-64;Novel-38*;Novel-74;Novel-112	GO:0032760	positive regulation of tumor necrosis factor production	4	biological process	0.0146
Novel miRNA	Novel-15;Novel-118;Novel-18b	GO:0019226	transmission of nerve impulse	3	biological process	0.0164
Novel miRNA	Novel-3;Novel-121*;Novel-98	GO:0007585	respiratory gaseous exchange	3	biological process	0.0285
Novel miRNA	Novel-19;Novel-118; Novel-70; Novel-122a; Novel-122b	GO:0072659	protein localization to plasma membrane	4	biological process	0.0347
Known miRNA	oar-miR-377-5p; oar-miR-376a-3p; oar-miR-410-3p; oar-miR-409-5p; oar-miR-218a; oar-miR-376a-5p; oar-miR-431; oar-miR-377-5p; oar-miR-381-5p; oar-miR-3956-3p	GO:0005737	cytoplasm	38	cellular component	0.0007
Known miRNA	oar-miR-376c-5p; oar-miR-23b; oar-miR-191	GO:0005844	polysome	3	cellular component	0.0226
Known miRNA	oar-miR-539-5p; oar-miR-377-3p; oar-miR-377-3p; oar-miR-377-5p	GO:0005623	cell	4	cellular component	0.0482
Conserved miRNA	miR-339a;miR-7c; miR-4454; miR-486*; miR-146b; miR-2285ak; miR-195*; miR-574c; miR-196a; miR-2284aa*;	GO:0005654	nucleoplasm	44	cellular component	0.0069
Conserved miRNA	miR-135a; miR-135b; miR-504; miR-196a; miR-196b; miR-7c; miR-767; miR-146b; miR-670; miR-574c; miR-339a; miR-874	GO:0016020	membrane	31	cellular component	0.0070
Conserved miRNA	miR-124; miR-215; miR-204a	GO:0031091	platelet alpha granule	3	cellular component	0.0289
Conserved miRNA	miR-135b; miR-208a; miR-1957; miR-146a; miR-2285ak; miR-7c; miR-331*; miR-383	GO:0043209	myelin sheath	8	cellular component	0.0292
Conserved miRNA	miR-486*; miR-1957; miR-208a; miR-345a; miR-2285ao; miR-2285ak; miR-339a; miR-670; miR-124; miR-144	GO:0005925	focal adhesion	13	cellular component	0.0458
Novel miRNA	Novel-5; Novel-121*; Novel-70; Novel-18b	GO:0098793	presynapse	4	cellular component	0.0118
Known miRNA	oar-miR-377-3p; oar-miR-412-5p; oar-miR-411b-5p; oar-miR-411b-5p; oar-miR-410-3p; oar-miR-376c-5p; oar-miR-539-5p; oar-miR-655-3p; oar-miR-409-5p	GO:0000166	nucleotide binding	9	molecular function	0.0040
Known miRNA	oar-miR-377-3p; oar-miR-376d; oar-miR-376c-5p	GO:0005262	calcium channel activity	3	molecular function	0.0271
Known miRNA	oar-miR-431; oar-miR-411b-5p; oar-miR-376a-5p	GO:0004722	protein serine/threonine phosphatase activity	3	molecular function	0.0376
Conserved miRNA	miR-2284d*-3; miR-670; miR-2284n*-4; miR-345a; miR-146a; miR-2285ak; miR-338*; miR-2284aa*, miR-2284d*-3, miR-2284aa*	GO:0000166	nucleotide binding	14	molecular function	0.0136
Conserved miRNA	miR-2284n*-4; miR-2284d*-3; miR-129; miR-486*; miR-146a; miR-2285ak; miR-542; miR-216; miR-2284aa*; miR-338*; miR-338*	GO:0044822	poly(A) RNA binding	29	molecular function	0.0388
Conserved miRNA	miR-628; miR-9*-1; miR-504; miR-2285ak; miR-208a; miR-574c; miR-216; miR-122; miR-339a; miR-29	GO:0005524	ATP binding	39	molecular function	0.0435

many miRNAs were in table, the first ten miRNAs were on the list

Table 11. Gene ontology analysis of differentially expressed genes in sheep skin between 75d and 85d of secondary HF morphogenesis.

miRNA category	miRNA	GO_id	GO_term	Count	Go category	P-Value
known miRNA	oar-miR-411b-5p;oar-miR-23b	GO:0048477	oogenesis	2	biological process	0.0441
Conserved miRNA	miR-2284n*-4;miR-32*;miR-133; miR-216; miR-574a; miR-338*;	GO:0000209	protein polyubiquitination	6	biological process	0.0147
Conserved miRNA	miR-2284aa*;miR-2284aa*;miR-216;miR-208a;miR-2284aa*;miR-2285av;miR-2285ak;miR-2285bb;miR-1;miR-206;	GO:0000398	mRNA splicing, via spliceosome	6	biological process	0.0301
Conserved miRNA	miR-2284aa*;miR-504*;miR-2284aa*; miR-2284d*-3;miR-2285af;miR-2285aq;miR-2285m*;miR-670; miR-196a;miR-196b; miR-133	GO:0006260	DNA replication	5	biological process	0.0332
Conserved miRNA	miR-15*;miR-504	GO:1,900,028	negative regulation of ruffle assembly	2	biological process	0.0384
Conserved miRNA	miR-32*;miR-133;miR-574a;miR-574d;miR-574f;miR-1;miR-338*	GO:0002223	stimulatory C-type lectin receptor signaling pathway	4	biological process	0.0490
Novel miRNA	Novel-18a;Novel-18b; Novel-7;Novel-20	GO:0002092	positive regulation of receptor internalization	3	biological process	0.0022
Novel miRNA	Novel-17;Novel-19; Novel-7	GO:0031122	cytoplasmic microtubule organization	3	biological process	0.0056
Conserved miRNA	miR-146b;miR-2285ak;miR-2285bb;miR-195*; miR-216;miR-1468;miR-146b;miR-628;miR-2284aa*;miR-2284aa*;miR-32*;miR-124;miR-2284n*-4;miR-216;miR-338*;miR-138*-1;miR-2284d*-3;miR-216;miR-196b;miR-208a;miR-2284aa*;miR-504*;miR-144;miR-208a;miR-2285av*-1;miR-196a;miR-196b; miR-1;	GO:0005654	nucleoplasm	32	cellular component	0.0067
Conserved miRNA	miR-195*;miR-216;miR-196a;miR-196b;	GO:0016581	NuRD complex	3	cellular component	0.0072
Conserved miRNA	miR-2284aa*;miR-338*;miR-2284d*-3;miR-2285af;miR-2285aq;miR-2284n*-4;miR-2284aa*;miR-2285av;miR-2285ak;miR-2285bb;miR-2284d*-3;miR-2285af;miR-2285aq;miR-2285m*; miR-15*	GO:0000166	nucleotide binding	8	molecular function	0.0144
Conserved miRNA	miR-628;miR-195*;miR-196a;miR-196b;	GO:0004407	histone deacetylase activity	3	molecular function	0.0332
Conserved miRNA	miR-628;miR-7b*-1;miR-504;miR-32*;miR-196a;miR-196b;miR-2285ak;miR-2285bb;miR-1468;miR-146b;miR-208a;miR-208b;miR-3432a;miR-34;miR-196b; miR-216;	GO:0005515	protein binding	76	molecular function	0.0356
Conserved miRNA	miR-125a;miR-133;miR-628;miR-7b*-1;miR-216;miR-208a;miR-2284n*-4;miR-2284d*-3;miR-2285af;miR-2285ak;miR-2285bb;miR-1;miR-206; miR-338*; miR-2284aa*;miR-2285av;miR-15*	GO:0044822	poly(A) RNA binding	15	molecular function	0.0375
Conserved miRNA	miR-138*-1;miR-2284d*-3;miR-2284aa*;miR-195*; miR-338*;miR-216;miR-133	GO:0000978	RNA polymerase II core promoter proximal region sequence-specific DNA binding	7	molecular function	0.0471
Conserved miRNA	miR-2284n*-4;miR-338;miR-504;miR-133;miR-196b*;miR-216	GO:0061630	ubiquitin protein ligase activity	5	molecular function	0.0492
Novel miRNA	Novel-64;Novel-51	GO:0005088	Ras guanyl-nucleotide exchange factor activity	2	molecular function	0.0325
Novel miRNA	Novel-19; Novel-7	GO:0008568	microtubule-severing ATPase activity	2	molecular function	0.0325
Novel miRNA	Novel-12;Novel-46;Novel-126;Novel-19;Novel-64;Novel-46;Novel-112;Novel-73; Novel-7;Novel-20	GO:0005524	ATP binding	10	molecular function	0.0479

many miRNAs were in table, the first ten miRNAs were on the list

Because different genes can coordinate their biological functions, in vivo KEGG pathway-based analysis can further elucidate gene function. In this study, we further analysed genes through the DAVID website (http://david.abcc.ncifcrf.gov/), and the results showed that 15 pathways are potentially involved in the corresponding regulation of miRNAs (Table 12). For example, mitogen-activated protein kinase (MAPK) and AMPK pathways regulate cell proliferation and differentiation, and the AMPK signalling pathway was the most significantly enriched (P-value = 0.0441). Differentially expressed genes in KEGG biological pathways were primarily involved in the MAPK signalling pathway, prolactin (PRL) signalling pathway, and in platelet activation, among other processes, which previous studies have shown to possibly be associated with biological pathways involved in HF development. Additionally, we conducted KEGG analysis of the target genes of differentially expressed miRNAs in sheep skin between 75 and 85 d of SF morphogenesis (Table 13), and a total of three enriched pathways were identified, including the PRL signalling pathway (P-value = 0.0145), platelet activation (P-value = 0.0292) and osteoclast differentiation (P-value = 0.0301). There were 2 common KEGG pathways in the 11 SF development periods, the PRL signalling pathway and platelet activation, and the most significant difference was in the PRL signalling pathway. These results showed that the PRL signalling pathway and platelet activation play important roles in SF development at 85 d. Five PRL signalling pathway genes (SOCS3, FOXO3, STAT1, PIK3R1, and CISH) were significantly enriched at 75 d compared with their expression at 85 d (corrected P-value<0.05). The genes STAT1 (STAT), CISH (Cis), and SOCS3 (SOCS) were significantly upregulated, and the gene FOXO3 was significantly downregulated (Figure 5). Six platelet activation genes (TLN2, FERMT3, RASGRP1, PIK3R1, PLA2G4D, and SYK) were significantly enriched at 75 d compared with their expression at 85 d (corrected P-value<0.05). The genes PLA2G4D (PLA), TLN2 (Talin), and FERMT3 (Kindlin3) were significantly upregulated, and the genes SYK, RASGRP1, and PIK3R1 (phosphatidylinositol-3-kinase, PI3K) were significantly downregulated (Figure 6). The target genes of the PRL signalling pathway and of platelet activation in sheep skin between 75 and 85 d of SF morphogenesis corresponded to the following miRNAs: oar-miR-133, Novel-17, miR-196b*, miR-2285 l, Novel-8*, oar-miR376c-3p, miR-504, Novel-64, miR-873, and Novel-20 (Table 14).

Table 12. KEGG Pathway analysis of target genes of differentially expressed miRNA in 11 periods development of sheep skin HF morphogenesis.

Pathway ID	Pathway	Gene name	miRNA	Gene counts	P-Value
ko04152	AMPK signaling pathway	PPP2R3A, PFKFB4, LEPR, PFKFB1, PPP2R5C, ADIPOR1, EEF2, FOXO3, PPARGC1A, CPT1A, AKT1, G6PC, GYS1, RAB14, EEF2K, PIK3R1, PPP2R2A	Novel-51；Novel-89；miR-196a；miR-196b；oar-miR-376a-5p；Novel-83；Novel-20；oar-miR-544-3p；miR-7c；Novel-17；oar-miR-409-5p；miR-105；Novel-13；miR-365；Novel-122a；Novel-122b；miR-144*；miR-7c；miR-2285ar；miR-2285l；oar-miR-376a-5p	17	0.0441
ko04921	Oxytocin signaling pathway	EGFR, CAMK1G, ADCY7, MAP2K2, PRKCG, EEF2, ITPR3, CACNG1, TRPM2, KCNJ9, EEF2K, RYR2, NFATC4, PPP3CA, PIK3R1, PLA2G4D	miR-486*；miR-383；oar-miR-377-3p；Novel-84；miR-4454；oar-miR-544-3p；miR-7c；Novel-111；Novel-121；Novel-3；Novel-5；miR-7c；miR-574c；Novel-74；miR-137；miR-2285ar；miR-2285l；miR-873；Novel-111；	16	0.0147
ko05223	Non-small cell lung cancer	AKT1, EGFR, FHIT, MAP2K2, PRKCG, FOXO3, STK4, PIK3R1	Novel-13；miR-486*；miR-365；Novel-84；miR-4454；Novel-17；miR-486；miR-7a*；miR-2285ar；miR-2285l；	8	0.0301
ko04917	Prolactin signaling pathway	AKT1, GCK, SOCS3, MAP2K2, GSK3B, FOXO3, STAT1, PIK3R1, CISH	Novel-13；miR-146a；oar-miR-133；Novel-84；oar-miR-26a；Novel-17；miR-196b*；miR-2285ar；miR-2285l；Novel-8*；	9	0.0332
ko04721	Synaptic vesicle cycle	TCIRG1, STXBP1, CLTC, UNC13C, RIMS1, ATP6V0D2, UNC13B, CACNA1A	miR-486*；Novel-70；oar-miR-3958-5p；oar-miR-218a； Novel-122a；Novel-122b；Novel-113； miR-29；Novel-18b	8	0.0384
ko04010	MAPK signaling pathway	AKT1, DUSP5, MAP2K2, RASGRP1, CACNG1, TRAF6, CACNA1A, PLA2G4D	Novel-13；Novel-74；Novel-84；Novel-64；Novel-121；Novel-112；Novel-18b；miR-873；Novel-111；	8	0.0490
ko04931	Insulin resistance	AKT1, G6PC, SOCS3, GSK3B, GYS1, TBC1D4, NFKBIA, PPARGC1A, PIK3R1, CPT1A, NR1H3	Novel-13；miR-365；oar-miR-133；oar-miR-26a；Novel-122a；Novel-122b；oar-miR-411b-3p；oar-miR-544-3p；oar-miR-409-5p；miR-2285ar；miR-2285l；miR-105；oar-miR-376c-5p	11	0.0022
ko04550	Signaling pathways regulating pluripotency of stem cells	AKT1, ACVR1B, TBX3, MAP2K2, WNT3A, GSK3B, DLX5, INHBC, ESX1, WNT6, ZFHX3, PIK3R1, AXIN1	Novel-13；Novel-82；oar-miR-3956-3p；Novel-84；Novel-118*；oar-miR-26a；oar-miR-377-3p；miR-184；miR-146b；Novel-121；oar-miR-381-3p；miR-2285ar；miR-2285l；miR-124；	13	0.0056
ko05213	Endometrial cancer	AKT1, EGFR, MAP2K2, GSK3B, FOXO3, PIK3R1, AXIN1	Novel-13, miR-486*, Novel-84, oar-miR-26a, Novel-17, miR-2285ar, miR-124	7	0.0067
ko04728	Dopaminergic synapse	AKT1, PPP2R3A, KCNJ9, DRD5, GSK3B, PPP2R5C, GRIN2A, PRKCG, PPP3CA, ITPR3, CACNA1A, PPP2R2A	Novel-13, Novel-51, Novel-5, Novel-15, oar-miR-26a, Novel-83, oar-miR-376d, miR-4454, miR-137, Novel-111, Novel-18b, oar-miR-376a-5p	12	0.0072
ko04611	Platelet activation	AKT1, FGG, ADCY7, TLN2, FERMT3, RASGRP1, RASGRP2, COL2A1, ITPR3, PIK3R1, PLA2G4D, SYK	Novel-13, miR-124, oar-miR-377-3p, oar-miR-376c-3p, miR-504, Novel-64, miR-574a, Novel-83, Novel-111, miR-2285ar, miR-873, Novel-20	12	0.0144
ko04920	Adipocytokine signaling pathway	AKT1, G6PC, SOCS3, LEPR, ADIPOR1, NFKBIA, PPARGC1A, CPT1A	Novel-13, miR-365, oar-miR-133, miR-196a, miR-196b, Novel-20, oar-miR-544-3p, oar-miR-409-5p miR-105	8	0.0332
ko05222	Small cell lung cancer	AKT1, FHIT, COL4A2, BIRC7, NFKBIA, TRAF6, BIRC3, PIK3R1, TRAF4	Novel-13, miR-365, Novel-121, Novel-13,oar-miR-544-3p,Novel-112, miR-204a, miR-2285ar, miR-2285l, Novel-98	9	0.0356
ko04973	Carbohydrate digestion and absorption	AKT1, G6PC, MGAM, SLC37A4, ATP1A1, PIK3R1	Novel-13, miR-365, miR-2285bd, miR-9*-1, miR-146a, miR-2285ar; miR-2285l	6	0.0375
ko05160	Hepatitis C	AKT1, EGFR, SOCS3, GSK3B, TICAM1, NFKBIA, PSME3, TRAF6, STAT1, PIK3R1, NR1H3, PPP2R2A	Novel-13, miR-486*, oar-miR-133, oar-miR-26a, Novel-38*, oar-miR-544-3p, miR-133, Novel-112, miR-196b*, miR-2285ar, miR-2285loar-miR-376c-5p, oar-miR-376a-5p	12	0.0471

Table 13. KEGG Pathway analysis of target genes of differentially expressed miRNA in sheep skin between 75d and 85d of secondary HF morphogenesis.

Pathway ID	Pathway	Gene counts	P-value
ko04917	Prolactin signaling pathway	5	0.0145
ko04611	Platelet activation	6	0.0292
ko04380	Osteoclast differentiation	6	0.0301

Table 14. Predicted target genes of prolactin signaling pathway and platelet activation in sheep skin between 75d and 85d of secondary HF morphogenesis.

Prolactin signaling pathway		Platelet activation
Target gene	microRNA	Target gene	microRNA
SOCS3	oar-miR-133	TLN2	oar-miR-376c-3p
FOXO3	Novel-17	FERMT3	miR-504
STAT1	miR-196b*	RASGRP1	Novel-64
PIK3R1	miR-2285l	PIK3R1	miR-2285l
CISH	Novel-8*	PLA2G4D	miR-873
		SYK	Novel-20

Figure 5. Differentially expressed genes between 75 d and 85 d of secondary HF morphogenesis involved in prolactin signalling pathway.

Note: The red-labelled genes were significantly upregulated at 75 d compared with expression at 85 d (P-value<0.05). The blue-labelled genes were significantly downregulated at 75 d compared with expression at 85 d (P-value<0.05). The green-labelled genes were expressed in sheep skin during secondary HF initiation, but there were no significant differences.

Figure 6. Differentially expressed genes between 75 d and 85 d of secondary HF morphogenesis involved in platelet activation.

Note: The red-labelled genes were significantly upregulated at 75 d compared with expression at 85 d (P-value<0.05). The blue-labelled genes were significantly downregulated at 75 d compared with expression at 85 d (P-value<0.05). The green-labelled genes were expressed in sheep skin during secondary HF initiation, but there were no significant differences.

Validation of small RNA sequencing by quantitative real-time PCR (qRT-PCR)

To validate the accuracy and reproducibility of the miRNAome results, we chose to present data for 4 of the 10 miRNAs that were significantly differentially expressed (miR-184 Novel-14 miR-206, miR-767) (Figure 7). The results of qRT-PCR analysis were partially consistent with the pattern shown by sequencing. miR-184 (1–2, 8–9), Novel-14 (5–6), miR-206 (1–2,2–3), miR-767 (3–4, 4–5, 8–9)did not show consistent expression between qRT-PCR and sequencing. The results for the rest of the miRNAs were consistent between qRT-PCR and sequencing analysis.

Figure 7. Validation of sequencing data by qRT-PCR.

The results showed that the qRT-PCR results were highly consistent with those of the sequencing data. Hence, RNA-seq provided reliable data. To clarify the miRNA abundance during the 11 developmental stages, qRT-PCR was used to analyse the expression pattern of oar-miR-133 (Figure 8).The expression level of this miRNA was highest at 85d (period 5), when it was significantly higher than duringother developmental periods (P-value<0.05), except for at 45d (period 1). The expression level of oar-miR-133 was significantly higher at 45 d (period 1), 55 d (period 2), 65 d (period 3), 75 d (period 4), 95 d (period 6), 105 d (period 7), and 125d (period 9) than at 115 d (period 8),135 d(period 10) and 145d (period 11; P-value<0.05). The results showed that oar-miR-133 plays an important role in SF morphogenesis at 85d.

Figure 8. Relative expression levels of oar-miR-133 in foetuses at different developmental stages.

Note: A different lowercase letter (a–c) above each bar indicates a significant difference among different developmental stages (P-value<0.05)

Discussion

Wool is an important animal product and raw material for the textile industry. The quality of wool is determined by theFD, crimp, yield, colour, and staple strength, among other properties.The quality and commercial value of animal hair depend on the structure and characteristics of HFs, and SFs determine the yield and FD of fine wool. The HF is one of the most unique and complex micro-organs in mammals, andis the only organ that exhibits cylindrical growth.

According to previous studies, after approximately 50 d of embryonic development, the skin of sheep begins to form PFs, and the first SFs begin to appear at up to E85 of development. At the foetal age of 105 d, secondary sub-follicles begin to form, and by birth, all HFs are formed, and PFs are mature. However, SFs do not fully matureuntil4–5 months after birth [23,24]. Additionally, studies have reported that PF development is initiated in sheep at E60, with follicle maturation occurring from E120, whereas SF differentiation occurs from E120 to the formation of cashmere after birth (NB). The foetal development period for PFs is E84 in Chinese super-fine-wool sheep (Gansu type), but in ourstudy,histomorphological observations of the foetal skin of Junken type sheep showed that the development of SFs occurred at E85, suggesting that the starting point of SF development differs among breeds.

miRNAs are involved in a variety of processes related to epidermal morphology and HF growth and development in animals. Some new methods have been used to studymiRNAs, especially high-throughput sequencing techniques. Zhang found that miRNAs are important in the function and evolution of skin based on a high-throughput sequencing techniques [25]. Li compared the miRNA transcriptomes of goat and sheep skin using Solexa sequencing to understand the development of skin and HFs and identified 1,910 known miRNAs and 2,261 novel mature miRNAs, of which 107 novel miRNAs and 1,246 known miRNAs were differentially expressed in the two libraries [26]. Emerging evidence indicates that regulatedmiRNA expression is essential for normal HF development [27–29]. Yi found that epidermal and HF miRNA expression is required for normal skin development in mice [30]. Mardaryev found that miRNA-31 expression in mouse HFs after the withdrawal period significantly increased with the administration of an antisense miR-31 inhibitory factor to the skin, which accelerated the development of HFs, cell differentiation and hair formation [31]. Moreover, investigations of HFs have provide invaluable insights into the possible causes of disorders related to human hair [32–34].

In this study,the miRNA expression profiles in the skin of sheep foetuses at11 developmental periods (45, 55, 65, 75, 85, 95, 105, 115, 125, 135 and 145d) were assessed by high-throughput sequencing and bioinformatics analysis.A total of 507 conserved miRNAs, 179 novel miRNAs and 152 known miRNAs were detected.72 conserved miRNAs, 44 novel miRNAs and 32 known miRNAs were significantly differentially expressed. The miRNAs showing differential expressionin the 11 developmental periods were clustered into 2 categories and 4 subcategories; one category included 45–95d of development, and the other included 105-145d of development.The first subcategory included 45d, 55d and 65d of development. The second subcategory included 75d, 85d and 95d of development. The third subcategory included 105d, 115d and 125d of development. The fourth subcategory included135d and 145d of development.The results of the cluster analysisare consistent with the findings for wool follicle development.Upon further analysis, oar-miR-10b showed the highest expression at 85 d of SF development.The formation and periodic development of mammalian HFs is under molecular regulation and requires the control of complex networks involving different signalling pathways. Liu et al found that differentially expressed miRNAs might regulate wool follicle development by targeting genes in many different pathways, such as the MAPK and Wnt pathways, as well as the pathways that regulate the actin cytoskeleton, focal adhesion, and tight junctions,during the anagen, catagen and telogenphasesin Tibetan sheep [35]. Zhang et al indicated that 7 miRNAs most likely target the genes in the Wnt/β-catenin, Shh/BMP and Notch pathways, which are important in feather morphogenesis. Based on our research, the AMPK and MAPK pathways, which are related to cell proliferation,differentiation and apoptosis,have an important link toHF development. KEGG pathway analysis showed that thetarget genes of the PRL signalling pathway and platelet activation were the most important among the differentially expressed miRNAs in SF morphogenesis.

Prolactin (PRL)has a direct effect on the division and proliferation of somatic cells [36]. Cassy showedthat melatonin acts indirectly on HF growth via the effects of PRL combined with its receptor (PRLR) [37]. PRL activates many signalling cascades, including the Jak2/Stat, which is the primary cascade, Src kinase, PI3K/AKT, and the MAPK pathways. Nixon et al. found that PRLR mRNA was most abundant in the dermal papilla, outer root sheath, hair germ, skin glands and epidermis. These cycle-related patterns of PRLR expression suggest dynamic regulation of PRLR by PRL that thereby modulates the hormonal responsiveness of seasonally growing HFs [38]. McCloghry found that PRL, but not melatonin, may be involved in the regulation of derived SF development [39]. Additionally, platelet activation and HFs are closely linked.Pietro found that platelet-rich plasma (PRP) has emerged as a new treatment modality in regenerative plastic surgery, and preliminary evidence suggests a beneficial role in hair regrowth [40]. Raquel analysed the adult HF dermal stem cell niche and suggested that platelet-derived growth factor signalling promotes self-renewal, which is essential for maintaining the dermal pool of HF stem cellsand, ultimately, their regenerative capacity within the HF [41]. In our study, the AMPK and MAPK pathways, which are related to cell proliferation, differentiation and apoptosis, [42,43] were found to have important links to HF development. The PRL signalling pathway and platelet activation have important rolesin SF morphogenesis. The genes SOCS3, FOXO3, STAT1, PIK3R1 and CISH were significantly enriched in the PRL signalling pathway (p < 0.05), and TLN2, FERMT3, RASGRP1, PIK3R1, PLA2G4D and SYK were significantly enriched in platelet activation (p < 0.05). Target genes in the PRL signalling pathway and of platelet activation corresponded to the following miRNAs: oar-miR-133, oar-miR-376c-3p, miR-196b*, miR-2285 l, miR-504, miR-873, Novel-8*, Novel-17, Novel-20 and Novel-64, which are candidate miRNAs for the initiation of SF development. In summary, our findings, combined with previous results, clearly indicated thatPRL and platelet activation have important roles in HF development. However, the mechanisms of action of signal transduction, apoptosis, nutrition and regulatory factors in HF development require further study.

Conclusions

E85 of foetaldevelopment is the starting point for SF development in Chinese Merino (Junken Type) sheep. High-throughput sequencingand bioinformatic analysis showed that a total of 507 conserved miRNAs, 179 novel miRNAs and 152 known miRNAs were detectedin 11skin expression profiles of sheep foetuses. A total of 37 differentially expressed miRNAs were detected in relation to SF morphogenesis. The PRL signalling pathway and platelet activation pathwaywereconfirmed to play important roles in SF development;oar-miR-133, oar-miR-376c-3p, miR-196b*, miR-2285 l, miR-504, miR-873, Novel-8*, Novel-17, Novel-20 and Novel-64 are potential candidate genes for further study of the molecular mechanisms underlying SF initiation.

Author Contributions

Hua Yang conceived and designed the study. Han-Yu-Lu Yang performed the experiments and wrote the manuscript. Hua Yang and Han-Yu-Lu Yang performed the statistical analysis and the data analysis. Xiao-Jun Liu,Min Shenand Jing–Quan Yangparticipated in performing the experiments. Guo-Qing Shi andYong-Lin Yangparticipated in the design of the study. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the the National Natural Science Foundation of China (31360541) [500000]; the Twelfth Five-Year-Plan of the Chinese 863 Program[2013AA102506] [1000000]; Guide Project of State Key Laboratory for Sheep Genetic Improvement and Healthy Production[SKLSGIHP2016A01] [350000].