Figures
Abstract
Cashmere goat skin contains two types of hair follicles (HF): primary hair follicles (PHF) and secondary hair follicles (SHF). Although multiple genetic determinants associated with HF formation have been identified, the molecules that determine the independent morphogenesis of HF in cashmere goats remain elusive. The growth and development of SHF directly influence the quantity and quality of cashmere production. Here, we report the transcriptome profiling analysis of nine skin samples from cashmere goats using 60- and 120-day-old embryos (E60 and E120, respectively), as well as newborns (NB), through RNA-sequencing (RNA-seq). HF morphological changes indicated that PHF were initiated at E60, with maturation from E120, while differentiation of SHF was identified at E120 until formation of cashmere occurred after birth (NB). The RNA-sequencing analysis generated over 20.6 million clean reads from each mRNA library. The number of differentially expressed genes (DEGs) in E60 vs. E120, E120 vs. NB, and E60 vs. NB were 1,024, 0 and 1,801, respectively, indicating that no significant differences were found at transcriptomic levels between E120 and NB. Key genes including B4GALT4, TNC, a-integrin, and FGFR1, were up-regulated and expressed in HF initiation from E60 to E120, while regulatory genes such as GPRC5D, PAD3, HOXC13, PRR9, VSIG8, LRRC15, LHX2, MSX-2, and FOXN1 were up-regulated and expressed in HF keratinisation and hair shaft differentiation from E120 and NB to E60. Several genes belonging to the KRT and KRTAP gene families were detected throughout the three HF developmental stages. The transcriptional trajectory analyses of all DEGs indicated that immune privilege, glycosaminoglycan biosynthesis, extracellular matrix receptor interaction, and growth factor receptors all played dominant roles in the epithelial-mesenchymal interface and HF formation. We found that the Wnt, transforming growth factor-beta/bone morphogenetic protein, and Notch family members played vital roles in HF differentiation and maturation. The DEGs we found could be attributed to the generation and development of HF, and thus will be critically important for improving the quantity and quality of fleece production in animals for fibres.
Citation: Gao Y, Wang X, Yan H, Zeng J, Ma S, Niu Y, et al. (2016) Comparative Transcriptome Analysis of Fetal Skin Reveals Key Genes Related to Hair Follicle Morphogenesis in Cashmere Goats. PLoS ONE 11(3): e0151118. https://doi.org/10.1371/journal.pone.0151118
Editor: Meijia Zhang, China Agricultural University, CHINA
Received: October 21, 2015; Accepted: February 22, 2016; Published: March 9, 2016
Copyright: © 2016 Gao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work was supported by grants from National Science Foundation Council (31372279, 31402038, and 31572369), Major Projects for New Varieties of Genetically Modified Organisms of China (2014ZX08008-002), Natural Science Foundation (2014JM3068) and International Sci-Tech Cooperation Program (2014KW14-01) of Shaanxi Province.
Competing interests: The authors have declared that no competing interests exist.
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.
(a) Schematic representation of the process of HF morphogenesis. Longitudinal sections (× 40) of fetal HF at (b) placode (E60), (c) differentiation stage (E120), and (d) maturation stage (NB). A high magnification view (× 160) of HF at (e) E120 and (f) NB. BL: Basal layer, DC: Dermal condensate, PHF: Primary hair follicle, SHF: Secondary hair follicle, DP: Dermal papilla, Bu: Bugle, SG: Sebaceous gland, HS: Hair shaft.
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.
Data from qPCR are shown as means ± standard error (SE) of three replicates. RPKM from RNA-seq are shown as means and SE of three replicates. The left side indicates the data from qPCR; the right side shows RPKM from RNA-seq.
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).
(a) E60: KLK11 was detected in the epidermis. The placode of HF shows little or no immunoreactivity at this stage (arrow). (b) E120: KLK11 was expressed strongly in the root sheaths and slightly in the bulb of HF (arrow). (c) NB: KLK11 was present in the hair matrix cells and root sheaths (arrow). (d) E60: VDR shows intense immunostaining in the epidermis. Very little VDR immunoreactivity was detected in the placode (arrow). (e) E120: VDR immunoreactivity is predominantly detected in the ORS keratinocytes and the bulb of HF (arrow). (f) NB: In fully formed HF, VDR staining was present in the root sheaths and bulb cells, where immunoreactivity was particularly enhanced in the ORS keratinocyte and matrix zone. Bar: 100 μm.
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).
The results are summarised in three main categories: biological process, cellular component and molecular function. The X-axis indicates the second level term of gene ontology; The Y-axis shows the percentage of genes.
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.
Profiles a-d: Each X-axis indicates the HF development state (E60, E120, and NB); The Y-axis shows expression changes. Trajectory cluster analyses (see Methods) show that gene expression did not steadily increase or decrease in the progression. Four profiles (profiles a-d) represent the major transcriptional trajectories.
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.
(a) The pathway and regulators during the development stages of HF initiation. (b) The pathway and regulators during development stages of HF cytodifferentiation and maturation. Arrows indicate either increased or decreased gene expression.
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.
https://doi.org/10.1371/journal.pone.0151118.s001
(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.
https://doi.org/10.1371/journal.pone.0151118.s002
(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.
https://doi.org/10.1371/journal.pone.0151118.s003
(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).
https://doi.org/10.1371/journal.pone.0151118.s004
(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.
https://doi.org/10.1371/journal.pone.0151118.s005
(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.
https://doi.org/10.1371/journal.pone.0151118.s006
(TIF)
S2 Table. Summary of read numbers based on the RNA-Seq data from cashmere goat HF development.
https://doi.org/10.1371/journal.pone.0151118.s008
(XLSX)
S3 Table. Gene with different expression in E60 vs. E120 and E60 vs. NB.
https://doi.org/10.1371/journal.pone.0151118.s009
(XLSX)
S4 Table. List of KEGG pathways for DEGs between different developmental stages.
https://doi.org/10.1371/journal.pone.0151118.s010
(XLSX)
S5 Table. List of genes with different expression in eight profiles.
https://doi.org/10.1371/journal.pone.0151118.s011
(XLSX)
Acknowledgments
We are grateful to Lei Qu and Jinwang Liu of Yulin University for the kind help and cooperation during the animal experiments.
Author Contributions
Conceived and designed the experiments: YG XW YJ YC. Performed the experiments: YG HY JZ SM YN. Analyzed the data: YG XW GZ. Contributed reagents/materials/analysis tools: YG HY JZ YN. Wrote the paper: YG XW YC.
References
- 1. Ansari-Renani HR, Ebadi Z, Moradi S, Baghershah HR, Ansari-Renani MY, Ameli SH. Determination of hair follicle characteristics, density and activity of Iranian cashmere goat breeds. Small Ruminant Research. 2011; 95(2–3):128–132.
- 2. Hardy MH. The secret life of the hair follicle. Trends in Genetics. 1992; 8(2):55–61. pmid:1566372
- 3. McDonald B, Hoey W, Hopkins P. Cyclical fleece growth in cashmere goats. Crop and Pasture Science. 1987; 38(3):597–609.
- 4. Galbraith H. Fundamental hair follicle biology and fine fibre production in animals. animal. 2010; 4(09):1490–1509.
- 5. Lee J, Tumbar T. Hairy tale of signaling in hair follicle development and cycling. Elsevier. 2012; 906–916.
- 6. Sennett R, Rendl M. Mesenchymal-epithelial interactions during hair follicle morphogenesis and cycling. Seminars in Cell & Developmental Biology. 2012; 23(8):917–927.
- 7. Yj Zhang, Yin J, Li Jq, Cq Li. Study on Hair Follicle Structure and Morphogenesis of the Inner Mongolian Arbas Cashmere Goat. Scientia Agricultura Sinica. 2007; 40(5):1017–1023.
- 8. Rendl M, Lewis L, Fuchs E. Molecular dissection of mesenchymal-epithelial interactions in the hair follicle. PLOS biology. 2005; 3(11):331.
- 9. Chang K-W, Huang NA, Liu IH, Wang YH, Wu P, Tseng YT, et al. Emergence of differentially regulated pathways associated with the development of regional specificity in chicken skin. BMC Genomics. 2015; 16(1):22.
- 10. Su R, Zhang WG, Sharma R, Chang ZL, Yin J, Li JQ. Characterization of BMP2 gene expression in embryonic and adult Inner Mongolia Cashmere goat (Capra hircus) hair follicles. Canadian Journal of Animal Science. 2009; 89(4):457–462.
- 11. Wu J-h, Zhang W-g, Li J-q, Yin J, Zhang Y-j. Hoxc13 Expression Pattern in Cashmere Goat Skin During Hair Follicle Development. Agricultural Sciences in China. 2009; 8(4):491–496.
- 12.
Krause K, Foitzik K. Biology of the hair follicle: the basics; 2006. WB Saunders. pp. 2–10.
- 13. Rogers GE. Hair follicle differentiation and regulation. International Journal of Developmental Biology. 2004; 48(2–3):163–170. pmid:15272381
- 14. Rishikaysh P, Dev K, Diaz D, Qureshi WMS, Filip S, Mokry J. Signaling Involved in Hair Follicle Morphogenesis and Development. International journal of molecular sciences. 2014; 15(1):1647–1670. pmid:24451143
- 15. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews Genetics. 2009; 10(1):57–63. pmid:19015660
- 16. Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq. Nature methods. 2011; 8(6):469–477. pmid:21623353
- 17. Geng R, Yuan C, Chen Y. Exploring Differentially Expressed Genes by RNA-Seq in Cashmere Goat (Capra hircus) Skin during Hair Follicle Development and Cycling. PLOS ONE. 2013; 8(4):e62704. pmid:23638136
- 18. Liu H, Wang T, Wang J, Quan F, Zhang Y. Characterization of Liaoning Cashmere Goat Transcriptome: Sequencing, De Novo Assembly, Functional Annotation and Comparative Analysis. PLOS ONE. 2013; 8(10):e77062. pmid:24130835
- 19. 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(1):30–38. pmid:23466980
- 20. Fan R, Xie J, Bai J, Wang H, Tian X, Bai R, et al. Skin transcriptome profiles associated with coat color in sheep. BMC genomics. 2013; 14(1):389.
- 21. Carter H, Clarke W. Hair follicle group and skin follicle population of Australian Merino sheep. Crop and Pasture Science. 1957; 8(1):91–108.
- 22. Auber L. VII.—The Anatomy of Follicles Producing Wool-Fibres, with special reference to Keratinization. Transactions of the Royal Society of Edinburgh. 1952; 62(01):191–254.
- 23. Nixon AJ. A method for determining the activity state of hair follicles. Biotechnic & histochemistry. 1993; 68(6):316–325.
- 24. Dong Y, Xie M, Jiang Y, Xiao N, Du X, Zhang W, et al. Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Capra hircus). Nature biotechnology. 2013; 31(2):135–141. pmid:23263233
- 25. Anders S, Huber W. Differential expression analysis for sequence count data. Genome biol. 2010; 11(10):R106. pmid:20979621
- 26. Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Annals of statistics. 2001:1165–1188.
- 27. Livak KJ, Schmittgen TD. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2− ΔΔCT Method. methods. 2001; 25(4):402–408. pmid:11846609
- 28. Langbein L, Schweizer J. Keratins of the human hair follicle. International review of cytology. 2005; 243:1–78. pmid:15797458
- 29. Giesen M, Gruedl S, Holtkoetter O, Fuhrmann G, Koerner A, Petersohn D. Ageing processes influence keratin and KAP expression in human hair follicles. Experimental dermatology. 2011; 20(9):759–761. pmid:21569108
- 30. Yu Z, Gordon SW, Nixon AJ, Bawden CS, Rogers MA, Wildermoth JE, et al. Expression patterns of keratin intermediate filament and keratin associated protein genes in wool follicles. Differentiation. 2009; 77(3):307–316. pmid:19272529
- 31. Rhee H, Polak L, Fuchs E. Lhx2 maintains stem cell character in hair follicles. Science. 2006; 312(5782):1946–1949. pmid:16809539
- 32. Laatsch CN, Durbin-Johnson BP, Rocke DM, Mukwana S, Newland AB, Flagler MJ, et al. Human hair shaft proteomic profiling: individual differences, site specificity and cuticle analysis. PeerJ. 2014; 2:506.
- 33. Hardy M, Vielkind U. Changing patterns of cell adhesion molecules during mouse pelage hair follicle development. Cells Tissues Organs. 1996; 157(3):169–182.
- 34. Kaplan ED, Holbrook KA. Dynamic expression patterns of tenascin, proteoglycans, and cell adhesion molecules during human hair follicle morphogenesis. Developmental dynamics. 1994; 199(2):141–155. pmid:7515726
- 35. Teti A. Regulation of cellular functions by extracellular matrix. Journal of the American Society of Nephrology. 1992; 2(10):S83. pmid:1318112
- 36. Powell BC, Nesci A, Rogers GE. Regulation of Keratin Gene Expression in Hair Follicle Differentiationa. Annals of the New York Academy of Sciences. 1991; 642(1):1–20.
- 37. Langbein L, Yoshida H, Praetzel-Wunder S, Parry DA, Schweizer J. The keratins of the human beard hair medulla: the riddle in the middle. Journal of Investigative Dermatology. 2009; 130(1):55–73.
- 38. Galbraith H. Protein and sulphur amino acid nutrition of hair fibre-producing Angora and Cashmere goats. Livestock Production Science. 2000; 64(1):81–93.
- 39. GRIMBLE RF. Sulphur Amino Acids. Nutrition and immune function. 2002; 1:133.
- 40. Iino M, Ehama R, Nakazawa Y, Iwabuchi T, Ogo M, Tajima M, et al. Adenosine stimulates fibroblast growth factor-7 gene expression via adenosine A2b receptor signaling in dermal papilla cells. Journal of Investigative Dermatology. 2007; 127(6):1318–1325. pmid:17301835
- 41. Pi L-Q, Jin X-H, Hwang ST, Lee W-S. Effects of calcitonin gene-related peptide on the immune privilege of human hair follicles. Neuropeptides. 2013; 47(1):51–57. pmid:22975462
- 42. Malgouries S, Thibaut S, Bernard B. Proteoglycan expression patterns in human hair follicle. British Journal of Dermatology. 2008; 158(2):234–342. pmid:18067481
- 43. Willen MD, Sorrell JM, Lekan CC, Davis BR, Caplan AI. Patterns of glycosaminoglycan/proteoglycan immunostaining in human skin during aging. Journal of investigative dermatology. 1991; 96(6):968–974. pmid:1710640
- 44. Sorrell J, Caterson B. Monoclonal antibodies specific for keratan sulfate detect epithelial-associated carbohydrates. Histochemistry. 1990; 94(3):269–275. pmid:1698187
- 45. Hubmacher D, Apte SS. The biology of the extracellular matrix: novel insights. Current opinion in rheumatology. 2013; 25(1):65. pmid:23143224
- 46. Akiyama M, Smith LT, Shimizu H. Changing patterns of localization of putative stem cells in developing human hair follicles. Journal of investigative dermatology. 2000; 114(2):321–327. pmid:10651993
- 47. Schneller M, Vuori K, Ruoslahti E. αvβ3 integrin associates with activated insulin and PDGFβ receptors and potentiates the biological activity of PDGF. The EMBO journal. 1997; 16(18):5600–5607. pmid:9312019
- 48. Lorenz K, Grashoff C, Torka R, Sakai T, Langbein L, Bloch W, et al. Integrin-linked kinase is required for epidermal and hair follicle morphogenesis. The Journal of cell biology. 2007; 177(3):501–513. pmid:17485490
- 49. 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(9):e76282. pmid:24069460
- 50. Inoue S, Nambu T, Shimomura T. The RAIG family member, GPRC5D, is associated with hard-keratinized structures. Journal of investigative dermatology. 2004; 122(3):565–573. pmid:15086536
- 51. Nachat R, Méchin M-C, Charveron M, Serre G, Constans J, Simon M. Peptidylarginine deiminase isoforms are differentially expressed in the anagen hair follicles and other human skin appendages. Journal of investigative dermatology. 2005; 125(1):34–41. pmid:15982300
- 52. Jave-Suarez LF, Winter H, Langbein L, Rogers MA, Schweizer J. HOXC13 is involved in the regulation of human hair keratin gene expression. Journal of Biological Chemistry. 2002; 277(5):3718–3726. pmid:11714694
- 53. Li J, Baxter RM, Weiner L, Goetinck PF, Calautti E, Brissette JL. Foxn1 promotes keratinocyte differentiation by regulating the activity of protein kinase C. Differentiation. 2007; 75(8):694–701. pmid:17459087
- 54. Ma L, Liu J, Wu T, Plikus M, Jiang T-X, Bi Q, et al. Cyclic alopecia'in Msx2 mutants: defects in hair cycling and hair shaft differentiation. Development. 2003; 130(2):379–389. pmid:12466204
- 55. Rendl M, Polak L, Fuchs E. BMP signaling in dermal papilla cells is required for their hair follicle-inductive properties. Genes & development. 2008; 22(4):543–557.
- 56. McDowall M, Edwards NM, Jahoda C, Hynd PI. The role of activins and follistatins in skin and hair follicle development and function. Cytokine & growth factor reviews. 2008; 19(5):415–426.
- 57. Watt FM, Estrach S, Ambler CA. Epidermal Notch signalling: differentiation, cancer and adhesion. Current opinion in cell biology. 2008; 20(2):171–179. pmid:18342499