Comparative Transcriptome Analysis of Fetal Skin Reveals Key Genes Related to Hair Follicle Morphogenesis in Cashmere Goats.

Introduction
Cashmere goats have double coats consisting of non-modulated fine inner hairs or cashmere fibres produced by secondary hair follicles (SHF) and guard hairs produced by primary hair follicles (PHF), which are invaginated into the basement membrane of the skin (epithelial and mesenchymal compartment) [1, 2]. Cashmere is a fine wool cashmere fibre (generally with diameter < 19 μm) that is used to produce soft luxurious apparel. The number and density of SHF, which affect the yield and diameter of the cashmere fibres, determines the value of the cashmere fleece [3]. Given that the generation of HF is established during early fetal life, and fibre characteristics are realised when the follicles are mature, examination of the processes and transcriptional regulatory mechanisms of the skin epithelium and skin appendage morphogenesis is therefore required to achieve maximum cashmere production [4].
HF morphogenesis is considered to occur in the major stages, for example during induction and initiation, differentiation and maturation [5]. The process of follicle morphogenesis has been studied extensively in murine models, but rarely in goats that produce fibres [5, 6]. A study on fetal HF morphogenesis of the Inner Mongolia Cashmere goats demonstrated that the hair placodes are formed at 55–65 days gestation (~E60), the SHF undergo rapid cytodifferentiation at 105–125 days gestation (~E120), and all PHF and some SHF mature at 135 days gestation [7]. HF morphogenesis is therefore a continuum process between 55 and 135 days of fetal life, which can represent the initiation, differentiation and maturation stages of HF.
The regulation of HF morphogenesis involves a series of complex molecular intercommunications between the single-layered epithelium and dermal cell condensate in skin. The epithelial-mesenchymal interface (EMI) during the organogenesis of HF is presumed to influence cell-substrate interactions [8, 9]. However, fewer studies have been carried out in animals whose skin is comprised of two types of HF, such as cashmere goats [10, 11]. The time point of fetal HF morphogenesis and related transcriptional gene expression in cashmere goats remain to be elucidated. A better understanding of the biological characteristics and regulation of HF morphogenesis may provide approaches to enhance the formation of fleece, whereby proper development is critically important for achieving maximum fleece production.
The structure of the HF in mammals is complex [12]. HF development takes place during fetal skin development and is modulated by extra-follicular macro-environmental factors [12–14]. Considering the complexity of the HF, studies of fetal skin have been valuable for fully identifying DEGs that appear to be developmentally regulated. RNA-seq is an unbiased technology approach for data collection [15, 16]. Previous findings on adult cashmere goats and sheep skin transcriptome analyses identified that many genes and pathways may be important for the regulation of HF cycling and coat colour [17–20]. These studies provided some insights into the genes that play versatile roles in the extra-follicular macroenvironment in goats and sheep. However, the extra-follicular macroenvironmental signalling that regulated fetal HF morphogenesis in cashmere goats is still not understood. Herein, a skin transcriptome analysis at three specific developmental stages, 60- and 120-day-old embryos and newborn (E60, E120 and NB, respectively), was conducted to examine the associated transcriptional genes governing the relevant processes. A number of stage-specific DEGs revealed here represent an important resource for identifying the molecular basis of fetal during skin and HF development in cashmere goats.

Materials and Methods

Animals and sample collection
Based on previous artificial insemination records (semen was collected from different rams), nine pregnant Shaanbei White cashmere goats (two years old, weighing 30–40 kg) at the same stage of pregnancy were selected from a breeding farm in Yulin, in the Shaanxi Province of China. All animals underwent identical feeding practices in accordance with the goat farm instructions. Six fetuses were delivered from six different females by caesarean at E60 and E120; E60 represented the initiation stage, and E120 represented the development stage. Three fetuses were also killed within two hours of birth, as a NB, from three different ewes; NB represented the PHF maturation stage. Each time point had three replicates. Skin samples from fetuses were collected from the right mid-side of the fetus, rinsed in ice-cold DEPC-treated water and cut into small pieces. Every skin sample was divided into two parts; one was stained and one was immediately frozen in liquid nitrogen for RNA-seq and qPCR analyses. All experimental procedures and the study design were carried out in accordance with the Care and Use of Laboratory Animals (Ministry of Science and Technology of China, 2006), and were approved by the Experimental Animal Manage Committee of Northwest A&F University under contract (NWAFU-31402038).

Skin staining
Skin tissues were prepared for histological sectioning, using the procedures Carter and Clarke (1957) [21] described. Skin samples from the different fetuses were placed in centrifuge tubes containing 4% paraformaldehyde solution (made with 0.1 M sodium phosphate buffer, pH = 7.4). Samples were then placed in small individual baskets and dehydrated through a series of graded ethanol. Processed skin samples were then embedded in paraffin, and serial vertical sections of skin were cut at a thickness of 5 μm. Sections were stained using the special tetrachrome stain ‘sacpic’ [22]. The sacpic method was previously used to determine the activity state of HF during postnatal periods [23].

Total RNA isolation, library construction and sequencing
Total RNA was isolated, using Trizol Reagent (Invitrogen) according to the manufacturer’s instructions, from each skin sample after grinding the sample in liquid nitrogen. The quality and concentration of the total RNA were determined using an Agilent 2100 Bioanalyzer (Agilent). RNA samples were stored at -80°C for later library construction and sequencing.
Nine RNA libraries for each skin sample were constructed, representing samples from the three time points during HF development. The libraries were as follows: E60_1, E60_2 and E60_3 as replicate libraries for the E60 experimental group, E120_1, E120_2 and E120_3 for the E120 experimental group, and NB_1, NB_2 and NB_3 for the NB experimental group. Oligo (dTs) were used to isolate poly (A) mRNA. The mRNA was fragmented and reverse transcribed using random primers. Second-strand cDNAs were synthesised using RNase H and DNA polymerase I. The double-strand cDNAs were then purified using the QiaQuick PCR extraction kit. The required fragments were purified via agarose gel electrophoresis and were enriched through PCR amplification. Finally, the amplified fragments were sequenced using Illumina HiSeq™ 2000 (GeneDenovo Co., Guangzhou, China) according to the manufacturer’s specifications. The raw sequencing data were submitted to the NCBI SRA database under access number SRP059481.

Mapping reads to the reference genome
The original sequencing-received image data were transferred into sequence data via base calling, which is defined as raw data or raw reads stored in the fastq format. Raw reads of all nine samples were pre-processed through the removal of containing adaptors-reads with more than 5% unknown nucleotides. Low-quality reads were also removed, in which the percentage of low-quality bases of quality value ≤ 5 was more than 50% in a read. The clean reads of each stage were aligned to the goat genome assembly CHIR_1.0 [24] using SOAP aligner/soap2. Mismatches of no more than two bases were allowed in the alignment, and uniquely mapped reads were obtained.

Expression annotation
For gene expression analysis, the number of unique-match reads was calculated and normalised to RPKM (reads per kb per million reads). Expression levels of each gene between two groups were compared to give an expression difference using DESeq, as Abders and Huber [25] described. The P value corresponded to differential gene expression at statistically significant levels [26]. FDR (False Discovery Rate) was used to determine the P value threshold. DEGs were defined as FDR ≤ 0.05 and absolute value of log2Ratio ≥ 1.
Gene expression data were normalised to 0, log2(E120/E60) and log2(NB/E60). The union set of all DEGs was used for further expression pattern analysis. STEM (Short Time-series Expression Miner, v1.3.8) was used to profile the gene set into eight expression patterns. Cluster profiles with a q value ≤ 0.05 were considered significantly expressed. GO (gene ontology) annotation was analysed using Blast2GO software (version 2.3.5) (https://www.blast2go.com/). Functional classification of the DEGs was performed using WEGO software. The KEGG (kyoto encyclopaedia of genes and genomes) pathway annotation was carried out using Blastall software against the KEGG database (http://www.kegg.jp/).

Confirmation of RNA-seq results with qPCR
First-strand cDNA was generated from 1 μg total RNA isolated from the skin samples using the Superscript™ First-strand Synthesis System (Invitrogen). To confirm the transcriptomic analysis results, eight genes were chosen consistent with those derived from sequencing, and were subjected to qPCR (quantitative reverse transcription PCR), including DNER (delta notch-like EGF repeat containing), PRSS35 (protease serine 35), HOXC13 (homebox C13), TRPV3 (transient receptor potential cation channel subfamily V member 3), BREVICAN, MYOG (myogenin), cystatin-M, and cytochrome P450 4X1-like. ACTB (beta-actin) was chosen as an internal reference gene since it had an equal RPKM value among the three stages. The qPCR was carried out on a CFX Connect Real-time PCR Detection System (Bio-Rad) using SYBR Premix Ex Taq (TaKaRa). The qPCR was run as follows: 50°C for 2 min and 95°C for 10 min, followed by 45 cycles of 95°C for 15 s, 60°C for 1 min, and 72°C for 45 s. Each qPCR analysis was performed in triplicate. Relative gene expression levels were calculated using the 2-△△Ct method [27]. The primers used for qPCR are listed in S1 Table.

Immunohistochemistry
Immunohistochemical staining was performed using a rabbit super-sensitive two-step immunohistochemical detection kit (PV-9001; Zhongshan Goldenbridge Biochemistry, Co., Ltd., Beijing, China). According to the manufacturer’s instructions, paraffin-embedded blocks were sectioned at 5 μm, and then were deparaffinized and rehydrated. After a boiling pre-treatment in the citrate buffer (pH 6.0) for antigen retrieval, slides were immersed in 3% hydrogen peroxide for 10 min to remove the endogenous peroxidase activity, blocked with blocking reagent in normal goat serum (ZLI-9020, ZSGB-BIO, Beijing, China) for 2 hours. They were then incubated at 4°C overnight in a humidified chamber with the following primary antibodies: VDR (vitamin D receptor) (rabbit, 1:300; Boster, Wuhan, Hubei, China), KLK11 (kallikrein 11) (rabbit, 1:100; Boster, Wuhan, Hubei, China). After being washed three times with PBS for 5 min each time, the samples were treated with a secondary antibody (PV-9002; Zhongshan Goldenbridge Biotechnology Co., Ltd., Beijing, China) at 37°C for 20 min. Then, the samples were washed three times with PBS for 5 min each time, and were treated with DBA for 1 min, counterstained with haematoxylin and countered under light microscopy.

Results

Establishment of the stages of HF morphogenesis
The time point and characteristics of HF morphogenesis during embryogenesis in cashmere goats has been reported to be between 55 and 135 days of fetal life, and this process is divided into three specific stages: initiation, differentiation and maturation [7]. Our histological study demonstrated that mesenchymal cells began to accumulate at E60 and the basal cells formed a visible hair germ (Fig 1a and 1b). At E120, PHF down-growth reached the subcutis and formed hair shafts at the upper end. A large number of SHF branched out around the PHF and a hair road was formed. A large number of SHF underwent rapid cytodifferentiation (Fig 1a, 1c and 1e). At NB, the PHF acquired their maximal length and prominent hair/cashmere shafts emerged through the epidermis. The PHF and a small amount of SHF were then matured (Fig 1a, 1d and 1f). In addition, there are no clear morphological structure differences of HF at E120 and NB. Our results are consistent with previous findings showing that the development of SHF demonstrated a hysteresis effect [7]. The process of HF development was found to have overlapping morphological states: HF initiation, differentiation of PHF and SHF, and the maturation of PHF.

Identification of expressed transcripts in the skin transcriptome
To quantify the gene expression patterns of fetal skin samples, we constructed nine cDNA libraries and then subjected them to deep sequencing using Illumina HiSeq2000. From each library, we obtained over 21.1 million raw reads, which was sufficient for a quantitative analysis of the gene expression patterns. After filtering the adaptor sequences, regions containing N sequences and low quality sequences, over 20.6 million clean reads were generated in each library. The percentage of clean reads among raw tags in each library ranged from 97.52% to 97.86% (S1 Fig), indicating that high quality RNA-seq data were obtained for further analysis. The clean transcripts obtained were then used for further analysis. An average of 15.3 million reads per sample was mapped to the goat genome (range: 12.8–19.5 million). Of the total reads, the rate of match reads was more than 62%, and the remaining reads were unmatched (S2 Table).
Gene coverage is the percentage of a gene covered by reads. The similarity distribution showed a comparable pattern with approximately 40% of sequences having a similarity of 80% from three biological replicates, suggesting good reproducibility of this method (S2 Fig). A correlation analysis of two parallel samples provided an evaluation of the reliability of experimental results and rationality of sampling. In our study, the scatter plot showed that the correlation values of two biological replicates at each stage were up to 0.92 (R2 ≥ 0.92) based on the RPKM values, suggesting sufficient reproducibility and rationality of sampling (S3 Fig).

Differentially expressed genes at E60, E120 and NB
To better examine the biological mechanism of HF morphogenesis, it is important to identify the DEGs at each stage. Mean RPKM values were generated out of three biological replicates of each experimental group. Expression levels of a distinct gene from two groups were compared to give an expression difference using the DESeq [25]. The number of DEGs in E60 vs. E120, E120 vs. NB, and E60 vs. NB was 1,024, 0 and 1,801, respectively, for transcripts detected with |log2
Ratio| ≥ 1.0 and q value ≤ 0.05. Furthermore, no DEGs in E120 vs. NB were identified, suggesting that the E120 and NB have similar transcript profiles. A total of 2,059 DEGs were found in E60 vs. E120 and E60 vs. NB (S3 Table).
Our analysis identified a set of genes belonging to keratin family member encoding genes (KRT) and keratin-associated protein encoding genes (KRTAP), which were markedly up-regulated in E120 vs. E60. Studies have shown that KRT and KRTAP are major structural proteins of the hair fibre and sheath, and their content is important for fleece quality [28, 29]. In addition, we compared our results with RNA transcriptomic data from 20–50 SHF or PHF in cashmere goats. Of these, 49 KRT genes and 30 KRTAP genes were annotated in the goat genome [24]. We found that more than half of KRT and KRTAP (21/31) were detectable in the HF of cashmere goats (Table 1). Most of the KRT and KRTAP genes were evolutionarily conserved, however the expression patterns in HF present some differences among different species, such as between humans and sheep, due to the distinctive features of hair and wool [30]. Thus, our results could be useful in finding the expression patterns of KRT and KRTAP in the HF of cashmere goats, whose composition and interactions could determine the fibre properties.
Moreover, the following genes were up-regulated: LRRC15 (leucine-rich repeat-containing protein 15), which may function in cell communication; LHX2 (LIM homeobox 2), which maintains stem cell characteristics in HF; VSIG8 (V-set and immunoglobulin domain-containing protein 8), which is enriched in the cuticle; and PRR9 (proline-rich protein 9), which serves as substrates for transglutaminases that are responsible for cross-linking (S3 Table) [31, 32]. This indicated that these DEGs are associated with structural protein biosynthesis during HF and skin development.

Validation of the sequencing data by quantitative PCR
To validate the accuracy and reproducibility of the transcriptome results, we selected eight genes for qPCR validation (Fig 2). Linear regression analysis of these gene expression ratios between RNA-seq and qPCR was highly correlated (r2 = 0.90) (S4 Fig). Hence, these qPCR results confirmed that our RNA-seq findings provided reliable data.

Immunostaining of KLK11 and VDR
The distribution of KLK11 and VDR was not detected in the placode at E60 (Fig 3a and 3d). The immunoreactivity increased as follicle development progressed. KLK11 immunoreactivity was localised to the nucleus of cells of the root sheaths and bulb of HF at E120 and NB (Fig 3b and 3c). VDR immunoreactivity was particularly enhanced in the ORS keratinocyte and matrix zone at E120 and NB (Fig 3e and 3f).

GO-term analysis of DEGs
To better understand the biological behaviour of HF and skin morphogenesis, 2,059 DEGs were categorised into three gene ontology categories: cellular component, biological process, and molecular function. DEGs between E60 vs. NB and E60 vs. E120 were categorised into 72 and 73 functional groups based on sequence homology. The top five functional categories of up-regulation DEGs in E60 vs. E120 included cell, single-organism process, cellular process, cell part and binding. The top five functional categories of up-regulation DEGs in E60 vs. NB included cell, cell part, single-organism process, cellular process and binding (Fig 4).
Under the cellular component category, a large number of up-regulation DEGs, as well as down-regulation DEGs, was categorised as cell part, cell and organelle in E60 vs. NB and in E60 vs. E120 (Fig 4). The cellular and developmental process in the biological process, binding and catalytic activity in molecular function were found to be greater in E60 vs. NB than in E60 vs. E120 (Fig 4). Further enrichment analysis was found to be related to the development process and signal transduction, and the detected DEGs were enriched in different terms related to HF and skin development. For example, the DEGs LHX2, BMP4 (bone morphogenetic protein 4), S100A1 (S100 calcium binding protein A1), ECM1 (extracellular matrix protein 1) and WIF1 (WNT inhibitory factor 1) were enriched in the transcriptional activator, secreted term, binding protein, glycoprotein and secreted protein, respectively.

KEGG pathway enrichment analysis of DEGs
KEGG analysis predicted that 55.5% (1,142/2,059) of DEGs were involved in 218 pathways. The top 20 KEGG pathways with the highest representation of DEGs are listed in Table 2. The maps with the highest DEGs representation were metabolic pathways (141 DEGs, 12.35%). “Pathways in cancer” (64 DEGs, 5.60%) was the first significantly highly enriched (q value ≤ 0.05). “Pathways in cancer” is a collection of many pathways and has an important functions in promoting cell proliferation and evading apoptosis. The cell adhesion molecules, axon guidance, pathogenic Escherichia coli infection, tight junction, phagosome, leukocyte transendothelial migration and ECM (extra-cellular matrix)-receptor interaction are also significantly enriched, respectively (S4 Table). From these pathways, information on communication between and within the epidermis and dermis, which regulates HF development, can be obtained. As an example, there are two pathways for epidermis and dermis communication: cell adhesion molecules and ECM-receptor interaction pathway. DEGs in the cell adhesion molecules pathway is to mould, relax or reinforce cell contacts in areas of increased HF morphogenetic activity, such as NCAM1 (neural cell adhesion molecule 1) and CDH1 (cadherin 1, type 1, E-cadherin epithelial) [33]. ECM-receptor interaction pathways are important to HF morphogenesis because they affect cell physiological activities, including proliferation and migration [34, 35]. These annotations provided a good platform for further research into understanding EMI in the HF and skin development process.

Expression profiling among the three developmental stages
To determine the primary gene expression trajectories, we further conducted a clustering analysis of all 2,059 DEGs. These genes could be clustered into eight profiles based on their expression modulation via STEM software, in which 1,975 were clustered into four profiles (P value ≤ 0.05), including two down-regulated patterns (profile1 and profile0) and two up-regulated patterns (profile6 and profile7) (S5 Fig). Profile6 and profile1 contained genes modulated after E120, and profile0 and profile7 contained genes positively or negatively modulated along the whole time course (S5 Fig). Profile1 and profile0 contained 574 and 449 DEGs, respectively, while profile6 and profile7 contained 812 and 140 DEGs (S5 Fig). The results offer new information related to further characterization of novel molecules associated with HF and skin development in cashmere goats, and their expression profiles were characterised in detail.
The most abundant group was profile6, with 1,203 genes up-regulated in E120 samples, compared to E60 samples, before reaching steady-state (Fig 5c). At E120, a number of new SHF grew rapidly as branches of mature PHF. These up-regulated genes may be associated with SHF cytodifferentiation and PHF maturation. KRT and KRTAP are expressed in the HF of sheep and humans, which are abundantly expressed during HF differentiation [28, 30, 36, 37]. A total of 30 up-regulated genes belonged to the KRT and KRTAP, such as KRT34, KRT33a, and KRT33b (S5 Table). These data were in agreement, showing that hair keratinocyte proliferation was active after E120, and these genes may be essential in determining the structure and quality of cashmere fibres. Amino acid metabolism was significantly important to protein synthesis. Sulphur-containing amino acids stimulate proliferation of hair forming cells in the HF, including methionine, cysteine and cystine [38]. In the cysteine and methionine metabolism, four genes were annotated to profile6, including AHCY (adenosylhomocysteinase) (S5 Table). AHCY is a ubiquitous enzyme that catalyses the hydrolysis of S-adenosylhomocysteine to adenosine and homocysteine [39]. Adenosine and homocysteine indicate an important branch point in the metabolism of methionine and cysteine, and adenosine could stimulate the expression of FGF7 (fibroblast growth factor 7) in dermal papilla (DP) cells [40]. Methionine is a nutritionally essential amino acid for cashmere fibre production, and cysteine is also required to produce maximum growth [38]. These DEGs in cysteine and methionine metabolism contributed to HF and skin protein metabolism in cashmere goats.
We also found that genes associated with the immune system and infectious diseases did not reach peak expression until HF morphogenesis was almost completed in profile7 (Fig 5d). Evidence in mice and humans indicates that the HF in the anagen phase is an immune-privileged organ, which expressed at very low levels of Class II MHC class antigens [41] (S5 Table). Expression of these genes provides strong evidence that immune privilege is crucial during PHF initiation.
In the metabolic pathways, six DEGs were clustered to profile0 or profile1 showing down-regulated trends. They encode enzymes that regulate GAG (glycosaminoglycan) biosynthesis-keratin sulfate, including B4GALT4 (UDP-Gal:betaGlcNAc beta 1,4-galactosyltransferase, polypeptide 4), B4GALT2, and B3GNT1 (UDP-GlcNAc:betaGal beta-1,3-N-acetylglucosaminyltransferase 1) (S5 Table). In rodent models, GAG distribution around the hair matrix is essential for hair morphogenesis [42, 43]. Keratin sulfate is expressed in the keratinocytes of epithelial tissue and the connective tissue sheath of adult HF [44]. All of the aforementioned evidence suggests an important role for keratin sulfate activation at the time of placode formation.
In the ECM-receptor interaction pathway, 16 out of 23 DEGs were clustered to profile0 or profile1, showing down-regulated trends (Fig 5a and 5b). ECM is composed of fibrous structural proteins, for example, collagens and laminin, and matricellular proteins, for example, thrombospondin and tenascins, which modulate EMI and HF morphogenesis [34, 45]. Many of these ECM genes were observed in our DEG datasets, including COL5A2 (collagen type V, alpha 2), THBS4 (thrombospondin 4) and TNC (tenascin C) (S5 Table). We assume that the enhanced expression of genes plays multiple roles in PHF initiation.

Discussion
The yield and diameter of fleece in certain mammal species are thought to be determined by the number and density of HF, which are established in the early stages of fetal life [3]. However, HF development is a complex process, and the processes and developmental mechanisms of fleece production at different stages in cashmere goats remain largely uncharacterised. Previous studies have confirmed that three specific developmental stages are typical (E60, E120 and NB) during the time interval spanning from initiation to maturation of HF using sacpic staining [7]. Therefore RNA-seq was conducted using skin tissues from both fetal and newborn goats. Fibre diameter is highly correlated with the size of the dermal papilla cells, whose origin can be traced to the earliest features of the developing HF.
HF initiation is thought to be partially influenced by interactions of ECM-receptor components [34]. ECM-receptors consist of integrins, proteoglycans and other transmembrane molecules that mediate the ECM interactions with cells [45]. Integrins have been shown to be expressed in mesenchymal aggregates of developing HF [46]. We found that three α-integrin genes were highly expressed at E60: α5, α9 and α11 (Table 3, S6 Fig). Integrins frequently act synergistically with other growth factor receptors [47, 48]. We also discovered that two corresponding growth factor receptors genes were highly expressed at E60, including FRS3 (fibroblast growth factor receptor substrate 3) and FGFR1 (fibroblast growth factor receptor 1) (Table 3, S6 Fig). These findings were in line with the DP and skin high-throughput transcriptome sequencing results obtained during mature HF cycling phases, indicating that these genes probably have a role in HF EMI, and they might be the mediators of HF initiation in cashmere goats [19, 49].
Some transcriptional regulator genes for HF protein synthesis were also remarkably up-regulated at E120 and NB, such as GPRC5D (G-protein coupled receptor family C group 5 member D) [50], PADI3 (protein-arginine deiminase type-3) [51], HOXC13 [52], FOXN1 (forkhead box protein N1) [53], and MSX2 (msh homeobox 2) [54] (Table 4, S6 Fig). This suggests that these genes take part in the keratinization process accompanied by basal regulation of the inner root sheath and hair shaft formation in cashmere goats. Here, we observed that Wnt-related genes were up-regulated at E120 and NB, including WNT10A (wingless-type MMTV integration site family, member 10A), WNT11, and β-catenin (Table 4, S6 Fig), implying that these genes probably play critical roles in HF cytodifferentiation and maturation. In addition, one other major signalling transduction pathway, TGF-beta (transforming growth factor-b)/BMP (bone morphogenetic protein), has recently been shown to play a central role in HF and skin development [14]. The Bmp4/Bmp2/Msx2/Foxn1 acidic hair keratin pathway is involved in the control of hair shaft growth and differentiation, and TGF-beta/BMP signals are necessary for regulating hair shaft growth and differentiation [5, 54, 55].
Expression of TGF-beta/BMP-related genes was also increased during the late stages of HF and skin morphogenesis, including BMP2, BMP4 and BMP8A. We also found that some BMP inhibitors were up-regulated in E120 skin, such as SMAD6 and SMAD7, which are important for antagonising TGF-beta/BMP activity and balance BMP inhibition [56] (Table 4, S6 Fig). Notch signalling controls DP signature molecules, which in return signal to the matrix cells to promote HF differentiation [5]. We noticed that up-regulated genes were related to Notch signalling, including JAG1 (jagged-1) (Table 4, S6 Fig). JAG1 is expressed in pre-cortex cells and the cuticle layer of the inner root sheath [57]. We believe that these annotated genes are involved in HF cytodifferentiation and maturation in cashmere goats.
Although Wnt and TGF-beta/BMP signalling pathways are all involved in HF formation and differentiation, we assume a mechanism by which HF development is regulated by key players within these pathways (Fig 6). For instance, BMP2 and BMP4 are down-regulated, and DKK1 is up-regulated during HF initiation.
In conclusion, our findings have greatly expanded the understanding of transcriptional responses within the three distinct developmental stages of HF. Novel expression patterns for thousands of genes were successfully established during HF morphogenesis. Furthermore, we hypothesise that some DEGs in the three signal transduction pathways (Wnt and TGF-beta/BMP) are independently involved in the cytodifferentiation and maturation of HF.

Supporting Information
S1 Fig
Classification of total raw reads at different developmental stages.
After filtering the adaptor sequences, regions containing N sequences and low quality sequences, the nine RNA-seq libraries generated over 20.6 million clean reads in each library. The percentage of clean reads among the raw reads reached 97.52% and 97.86% in each library.
(TIF)

S2 Fig
Percent of coverage representing the percentage of a gene covered by reads at each stage.
The distribution of distinct reads over different read abundance categories showed similar patterns for all nine RNA-seq libraries. The similarity distribution showed a comparable pattern with approximately 40% of the sequences having a similarity of 80% from the three biological replicates.
(TIF)

S3 Fig
Correlations between biological replicates and skin samples.
The x- and y-axis correspond to the RPKM value of each sample. The correlation coefficient (r2) between two individuals within each group was calculated based on the RPKM value of each individual. Correlation values of two biological replicates at each stage were up to 0.90.
(TIF)

S4 Fig
Coefficient analysis of fold change data between qPCR and RNA-seq.
Eight genes were selected for qPCR. Data indicating relative transcript level from qPCR and RPKMs from RNA-Seq are means of three replicates in each group. Scatterplots were generated by the log2expression ratios from RNA-seq (X-axis) and qPCR (Y-axis).
(TIF)

S5 Fig
Patterns of gene expression by STEM across three phases.
Each profile represents an expression pattern. Patterns with colours indicate genes were significantly enriched in this pattern, while blank ones represent non-significance. The number of genes belonging to each pattern is labeled above the profile.
(TIF)

S6 Fig
Hierarchical cluster analysis of gene expression profiles from nine skin samples with 50 DEGs.
Columns are clustered by libraries and rows are clustered by genes. Dendrogram height indicates distances between clusters in gene expression profiles. Orange indicates up-regulation and blue indicates down-regulation. There were clusters with relatively minor differences for E120 vs. NB. The bottom of each column indicates three replicates in each HF development stage, from left to right, starting with E60_1, E60_2, E60_3, E120_1, E120_2, E120_3, NB_1, NB_2, and NB_3.
(TIF)

S1 Table
List of primers for qPCR.
(XLSX)

S2 Table
Summary of read numbers based on the RNA-Seq data from cashmere goat HF development.
(XLSX)

S3 Table
Gene with different expression in E60 vs. E120 and E60 vs. NB.
(XLSX)

S4 Table
List of KEGG pathways for DEGs between different developmental stages.
(XLSX)

S5 Table
List of genes with different expression in eight profiles.
(XLSX)