Transcriptomic analysis reveals critical genes for the hair follicle of Inner Mongolia cashmere goat from catagen to telogen

There are two main types of hair follicle in Inner Mongolia Cashmere goats, the primary hair follicle (PHF) producing hair fibers and the secondary hair follicle (SHF) producing cashmere fibers. Of both fibers from cashmere-bearing goats in Aerbasi, Inner Mongolia, the timing of cyclical phases for the cashmere have been well clarified but hair fibers have been less noticeable. Herein, we evaluated transcriptome of PHF and SHF from the same three goats in Aerbasi at the catagen- and telogen phase of cashmere growth. We totally found 1977 DEGs between PHFs at the telogen and catagen phases of SHF, 1199 DEGs between telogen- and catagen SHF, 2629 DEGs between PHF at the catagen phase of SHF and catagen SHF, and 755 DEGs between PHF at the telogen phase of SHF and telogen SHF. By analyzing gene functions based on GO and KEGG database, we found that the DEGs have functions in muscle contraction and muscle filament sliding between catagen- and telogen SHF, indicating that arrector pilli muscles might play a role on the transition from catagen to telogen. Moreover, considering that the enriched GO and KEGG categories of the DEGs between PHF and SHF, we suggested that part of PHF might rest in their own anagen phase when SHF are at catagen, but PHF might enter into the telogen phase at SHF’s telogen. Finally, we highly recommended the several potential genes acting as the regulators of the transition between growth phases including IL17RB and eight members of ZNF. These results provide insight into molecular mechanisms on the transition of SHF from catagen to telogen together with PHF’s growth situation at SHF’s catagen and telogen in Inner Mongolia Cashmere goats.


Introduction
Hair follicles on animal skin execute many useful biology function except for producing hair, such as immune defense and grease secretion, etc. There are two main types of hair follicle in goats, the primary hair follicle (PHF) producing hair fibers and the secondary hair follicle (SHF) producing cashmere fibers. Cashmere is finer, stronger, lighter, softer, and more insulating than sheep wool, surely it looks as a luxury fiber [1]. As a major province of raw cashmere in China, Inner Mongolia bears three types of Cashmere goats-Erlangshan, Aerbasi and Alashan; these Cashmere goats are mainly distributed in the northwest of Inner Mongolia Plateau (Fig 1). They are superior in the production of cashmere with diameter less than 16 μm.
Hair follicles is a unique, highly regenerative organ with a normal developmental cycle, which is mainly determined by the interaction of epithelial cells and dermal papilla that are generally composed of dermal mesenchymal cells [2]. A growth cycle consists mainly of three distinct stages-anagen, catagen and telogen [2,3,4]. In Cashmere goat, SHF takes approximately one year to complete a growth cycle. For Inner Mongolia Cashmere goats, the anagen phase of SHF begins from April until November, the catagen phase is from December to January, and the telogen phase is from February to March [5]. HFs help cashmere to grow at the anagen phase, turn into apoptosis at the catagen phase and finally release cashmere at the telogen phase. When the next anagen phase comes, HFs launch a new regenerative cycle. To increase the harvest accompanying with the thinner cashmere, researchers try to find major genes and pathways, which may affect the of development and growth the cashmere by molecular biology methods.
RNA-seq is usual and convenient for researchers to study different animal phenotypes due to differentially expressed genes (DEGs), because we could get enough transcriptome data from small quantities of tissue samples [6,7]. In previous studies of hair follicle in Cashmere goats, the scientists have already found many important factors that may affect hair follicle cycle in Wnt signal transduction pathway, fibroblast growth factor (FGF) family, bone morphogenetic protein (BMP) family, Sonichedgehog (Shh) signal transduction pathway, transforming growth factor (TGF) family, Notch signal transduction pathway and so on [8,9,10]. Some genes play an activator role while other genes play roles in these pathways as the inhibitors. β-catenin has been confirmed as an important accelerator in Wnt signal transduction pathway, but more importantly, its activity is simultaneous with the apoptosis of hair follicle stem cells [11]. Maksim V. Plikus found that BMPs may be the long-sought inhibitors of hair growth postulated by classical experiments [12].
In this study, we used RNA-seq technology to find some important factors in SHF of Inner Mongolia Cashmere goats by comparing catagen SHF with telogen SHF. Furthermore, functional annotation analysis also helps us find some important genes that influence cashmere growth. By comparing the DEGs of PHF at the catagen and telogen phase of SHF, the growth condition of PHF became relatively clear. Our data might take insight into the molecular mechanism on the growth of SHF and PHF at the catagen-and telogen phase of SHF.

Ethics statement
In this study, hair follicles were collected in accordance with the International Guiding Principles for Biomedical Research Involving Animals and was approved by the Animal Ethics Committee of the Inner Mongolia Academy of Agriculture and Animal Husbandry Sciences that is responsible for Animal Care and Use in the Inner Mongolia Autonomous Region of China. In our study, no specific permissions were required for these activities and the animals did not involve endangered or protected species.

Hair follicle samples preparation for RNA-seq and qRT-PCR validation
Three 3-years-old female Aerbasi Inner Mongolia Cashmere goats from goat stud farm (Aerbasi White Cashmere Goat Breeding Farm, Erdos, Inner Mongolia) were used in this study for RNA-seq. At the dorsal side of goats, we collected PHF and SHF samples by pulled out hair and cashmere from the root of HFs in SHF's catagen (in mid-January) and its telogen (in mid-March), respectively; and then put them in liquid nitrogen as soon as possible and finally stored in -80˚C refrigerator for a long term. About 50 to 100 hair follicles of three cashmerebearing goats were collected as a sequencing sample. Six female goats, which come from the same group at same age as the RNA-seq goats, were used in qRT-PCR.

Construction of RNA library, sequencing and quality control
Total RNA was extracted and test the concentration. Each RNA sample needs over than 2 μg and RIN (the integrity of RNA) value larger than six. We mixed three cashmere goats total RNA and then constructed a library. In hence, four libraries were constructed. After high throughput sequencing, we used Fast-QC to measure the quality of sequence data, including base quality distribution, GC%, content of PCR duplication and frequency of kmer.

Expression level analysis
RPKM means the genes expression and calculating bases on the reads counts mapped to this gene and gene length [14]. The different expression genes between catagen and telogen was calculated by EBSeq (Log2FC>1 or Log2FC<-1, FDR<0.05) [15].

GO and KEGG pathway analysis
We analyzed genes function by Gene Ontology (GO), an international gene functional classification system, on three categories including biological process, molecular function and cellular component [4,16]. We used Fisher formula to exam the P-Value of each GO analysis and blasted result onto human and mouse reference genomes (P-Value<0.01 represent gene function is significance). Different expression genes (DEGs) were annotated by Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg) to identify the pathway enrichments [17,18,19].

qRT-PCR to validate RNA-seq data
To validate the veracity of RNA-seq data and genes expression level, we used quantitative real time PCR in this experiment. First, we extracted total RNAfrom catagen's-and telogen's SHF from six goats. Second, we synthesized cDNA from mRNA. Then, the primers we used was designed and synthesized depending on the mRNA sequences published on NCBI database. QRT-PCR was used SYBR Green II. In the qRT-PCR, β-actin was acted as internal reference and DEGs expression level was calculated by 2 -ΔΔct . The primer sequences information was shown in Table 1.
Reaction system of qRT-PCR was totally 20μl, including SYBR Green II 10μl, primer F 0.5μl, primer R 0.5μl, ROX 0.4μl, cDNA 2μl and H2O 6.6μl. The condition of the reaction was departed into amplification stage and dissociation stage. Amplification stage includes 95˚C 30s for 1 cycle, 95˚C 5s and annealing temperature (TM) 34s for 40 cycles, dissociation stage includes 95˚C 15s, 60˚C 60s and 95˚C 15s. Result were analyzed by SPSS 17.0.

Result
To gain insight into the growth of hair follicles at SHF's cyclic phases, we quantified the wholegenome transcriptomes of primary-or secondary hair follicles at the catagen and telogen phases of SHF in 3 female Inner Mongolia Cashmere goats. In total, we obtained over 25 million clean reads from each library after trimming, adaptor sequences, low quality reads and multiple mapped reads ( Table 2). Gene structure was showed in Fig 2. The clean reads were mapped to 22175 genes annotated in the goat reference genome (Capra Hircus 1.0). By comparing the RNA-seq data, we then identified DEGs in the following four groups including (i) PHF at the telogen and catagen phases of SHF, (ii) telogen-and catagen SHF, (iii) PHF at the catagen phase of SHF and catagen SHF, and (v) PHF at the telogen phase of SHF and telogen SHF (Table 3). Using the thresholds of false discovery rate (FDR) < 0.05 and difference ratio of RPKM (reads per kilobase of exon per million fragments mapped) > 2, we defined 1977 DEGs between PHF at the telogen and catagen phases of SHF, 1199 DEGs between telogen and catagen SHF, 2629 DEGs between PHF at the catagen phase of SHF and catagen SHF, and 755 DEGs between PHF at the telogen phase of SHF and telogen SHF. When focusing on the characterized transcripts, the number of DEGs among the four groups arrays in ascending order as follows: 1783 DEGs between PHF at the catagen phase of SHF and catagen SHF, 1092 between PHF at the telogen and catagen phases of SHF, 552 between telogen and catagen SHF, and 422 between PHF at the telogen phase of SHF and telogen SHF. It is noteworthy that, based on the number of DEGs, the largest distinction at transcription level should be between PHF at SHF's catagen phase and catagen SHF in sharp contrast to the smallest difference between PHF at SHF's telogen phase and telogen SHF. Interestingly, the transcriptomic  changes seem to be relatively small between telogen-and catagen SHF compared with those between PHF at the catagen-and telogen phase of SHF. Next, we performed a GO annotation analysis for these characterized DEGs in the four groups. The biology processes were considered to be significant with the thresholds of false discovery rate (FDR) < 0.05 as showed in Table 4. Hairs detach dermal papilla and stop growing since catagen, and then they become dull and lifeless at telogen because of no blood supply. For the only muscle in hair follicles, the arrector pilli muscles is not just an attacher to hair follicle-as had been reported earlier-but is an essential part of the hair follicle cycle by interacting with the follicle mesenchyme [20]. As expected, between catagen SHF and telogen SHF, the DEGs were related to muscle contraction (GO:0006936) and muscle filament sliding (GO:0030049) instead of growth and metabolism. Our results indicated that arrector pilli muscles, which attached to hair follicles, might have functions on helping cashmere stay in the hair follicles at catagen. Interestingly, the DEGs of PHF between the two phases of SHF were involved in those biology processes including translation (GO:0006412, 0006413, 0006414 and 0006415), RNA metabolism (GO:0016070 and 0016071), protein metabolism (GO:0044267), and viral growth (GO:0019083, 0019058 and 0016032). There is not any enriched GO category of the DEGs between PHF and SHF at the telogen phase of SHF, which is satisfied with the threshold at the telogen phase of SHF. In contrast, we found 22 enriched GO categories for the DEGs between the PHF at the catagen of SHF and catagen SHF, such as extracellular matrix organization (GO:0030198), muscle contraction (GO:0030049 and 0006936), translation (GO:0006413, 0006414, and 0006415), and metabolism (GO:0044281, 0006749, 0005975, and 0006805). All in all, these results might hint that the catagen phase of PHF might overlap with the catagen phase of SHF.
Besides the GO analysis, we also enriched the KEGG pathways for DEGs with FDR < 0.05 in order to collect the molecular interaction, reaction and relation among DEGs. As same as the enriched biological processes in the GO analysis, the pathways for the DEGs between catagen SHF and telogen SHF are also associated with the function of muscle (Table 5). Furthermore, we also observed that the most pathways were observed for the DEGs between the PHF at the catagen of SHF and catagen SHF, including extracellular matrix interaction (PATH:04512 and 04510), metabolism (PATH:04974, 0006749, 00480, 00980 and 05240) and ribosome (PATH:03010). Against that several pathways for the comparisons as above described, the only enriched pathway for the DEGs of PHF between the two phases of SHF was ribosome (PATH:03010); besides, none can be enriched for the DEGs between PHF and SHF at the telogen phase of SHF.
According to GO terms, we made an examination using qRT-PCR for the DEGs, such as IL17RB (Interleukin 17 receptor B), HPS6(biogenesis of lysosomal organelles complex 2 subunit 3) and ALPL(alkaline phosphatase), which are functionally related to hair follicle growth Fig 3). The result showed that these three genes were upregulated in catagen SHF compared with telogen SHF, which is consistent with the result of RNA-seq. We also found that the DEGs between catagen SHF and telogen included eight members of zinc finger protein family (Table 6), of which seven members were upregulated except for ZNF347.  Critical genes for secondary hair follicle of cashmere goat

Discussion
In this study, we chose hair follicle as an alternative to skin in order to reveal the molecular mechanism of the growth cycle of Inner Mongolia Cashmere goat hair follicle. Compared with skin, hair follicle is beneficial to remove the signals of hair stem cells in skin and to eliminate some unexpected influence. To date, we are acquainted with the growth law of SHF in Inner Mongolia Cashmere goat but are completely at a loss for the time-based growth pattern of PHF yet. This study might be helpful to find the relationship between these two hair follicles by collecting PHF at the catagen/telogen phase of SHF. We uncovered that the distinction of gene expression are mainly associated with muscle movement between catagen SHF and telogen one. However, the growth cycle of PHF is not as same as SHF's; at SHF's catagen phase, PHF looks like stay in an active state because many metabolism-related genes were enriched (S1 Table). But, interestingly, PHF might be of the same inactive state as SHF at SHF' telogen phase because no one GO category was identified for DEGs between PHF and SHF. In addition, we analyzed the expression pattern of several critical genes in SHF at the catagen-and telogen. We found that a critical gene IL17RB was upregulated in catagen SHF compared with telogen SHF. IL17 superfamily is a relatively new family of cytokines, which consists of 6 ligands and can bind 5 receptor subtypes [21]. According to recent researches and our transcriptome data, we predicted that IL17RB might play a role in the development of SHF because it acts as a mediate signaling in many other physiological functions [22,23]. A branch of IL-17 pathway, IL17E pathway, targets epithelial cell. The compound of IL17RB and IL17RA worked as membrane receptor signals of IL17E which would activate TRADD, FADD and Casp in order. Finally, this pathway may indirectly affect the apoptosis of epithelial cells (Fig  4). It might indicate that IL17RB plays a positive role in cell apoptosis, and then launches the transition from catagen to telogen. Also, two IL17RB can also form a compound which acted as a membrane receptor of IL17B in IL17RB-IL17B pathway. Vahideh Alinejad et al. found that the IL17RB-IL17B pathway affected breast cancer together with the upstream and downstream cytokines [24]. They summarized IL17 family includes six protein members-IL17A, IL17B, IL17C, IL17D, IL17E and IL17F, of which IL17B and its receptor(IL17RB) are least known genes [25,26]. Then, we identified that the eight members of zinc finger protein family were the DEGs between catagen SHF and telogen one (Table 6), of which seven members were upregulated except for ZNF347. Zinc finger proteins are important structural elements in many nucleic acid binding proteins [27]. The researchers further found that zinc finger proteins might affect Critical genes for secondary hair follicle of cashmere goat the growth cycle of hair follicle [28,29]. A zinc finger transcription factor Trps1 mainly expressed at the nuclei of mesenchymal cells during hair follicle morphogenesis; its expression level increased in dermis but decreased in epidermis during early skin morphogenesis [28]. In 2013, Zfp157, a member of KRAB zinc finger protein family, expressed in both epithelial cells of ducts and sebaceous glands of hair follicles in mice [29]. Notably, ZNF667 can inhibit the expression and promoter activity of the rat proapoptotic gene Bax [30]. ZNF667 was a DEG from catagen to telogen of both SHF and PHF in our study, this may indicate that ZNF667 might be an important impact factor for the transition of SHF and PHF from catagen to telogen.
In our previous study, we found that STC2 was downregulated from anagen to catagen [31]. Teng Xu et al. identified approximate 7,000 transcripts that were differentially expressed between anagen and telogen of Inner Mongolia Cashmere goat. In addition, the genes were mainly enriched in ECM receptor interaction, focal adhesion and gap junction from the KEGG pathway database [32]. Rongqing Geng et al. totally identified 1,332 DEGs in Shaanbei Cashmere goat among anagen, catagen and telogen and these genes were mainly playing important roles in Wnt, Shh, TGF-β, and Notch signaling pathways [8]. Xiao-yang Ji et al. found that ubiquitin-mediated proteolysis pathway is a prominent signaling pathway that can distinguish SHF from PHF in cashmere goats [33]. Jianping Li et al. found ten microRNAs have influence in hair follicle growth in Liaoning Cashmere goat and fine-wool sheep telogen skins [34]. The researchers also kept on excavating inhibitors that may affect the growth of hair follicle from model animals to cashmere goats. Plikus et al. indicated BMPs might be a potential signal that played as an inhibitor in the propagation of hair stem cells in mice [12]. Wen L. Bai et al. found that BMP4 is generally upregulated from anagen to telogen in Liaoning Cashmere goat skin tissue, a higher BMP4 methylation level in skin coincides with a lower expression of BMP4 mRNA [35]. They also found seven downregulated miRNAs could joined some signal pathways through their target genes directly or indirectly from anagen to telogen [36]. Mathieu Castela et al. suggested IGF1R would together with BMP4 to control hair follicle goes into anagen early and into catagen late using K15-IGF1RKO mice [37]. Xianghui Ma et al. found MSI2, a RNA-binding protein of Musashi, would keep hair follicle resting in process of telogen-to-anagen transition in mice [38]. In recent study, Guangxian Zhou et al. revealed ncRNA might play a regulatory role in skin of cashmere goat through whole transcriptome [39]. Our study might improve the understanding on the molecular mechanism of the growth cycle of SHF and PHF in Cashmere goats.

Conclusion
In this study, we aim to find out about some interesting factors and signaling pathways, using transcriptome sequencing, that have functions on the growth of two kinds of hair in Inner Mongolia Cashmere goats at the catagen-and telogen phase of SHF. We totally found 1977 DEGs between PHF at the telogen and catagen phases of SHF, 1199 DEGs between telogen and catagen SHF, 2629 DEGs between catagen SHF and PHF at the catagen phase of SHF, and 755 DEGs between telogen SHF and PHF at the telogen phase of SHF. Between catagen-and telogen SHF, the DEGs were related to muscle contraction (GO:0006936) and muscle filament sliding (GO:0030049). Moreover, according to GO category of the DEGs between PHF and SHF, part of PHF might be in anagen phase at SHF's catagen, but PHF might enter into telogen phase at SHF's telogen. Among DEGs, IL17RB and eight members of ZNF might play important roles in the growth at the transition of hair follicles from catagen to telogen. These results would benefit the studies of the growth of two kinds of hair follicles in Inner Mongolia Cashmere goats.
Supporting information S1