Identification of potential key genes and pathways associated with the Pashmina fiber initiation using RNA-Seq and integrated bioinformatics analysis

Pashmina goat (Capra hircus) is an economically important livestock species, which habitats the cold arid desert of the Ladakh region (India), and produces a princely animal fiber called Pashmina. The Pashmina goat has a double coat fleece as an adaptation to the very harsh cold winters the outer long coarse hair (guard hair) produced from primary hair follicles and the inner fine Pashmina fiber produced from secondary hair follicles. Pashmina fiber undergoes a circannual and synchronized growth cycle. In the present study, we analyzed transcriptome profiles from 10 different Pashmina goats during anagen and telogen to delineate genes and signaling pathways regulating active (anagen) and regressive (telogen) phases of the follicle growth. During anagen, 150 genes were expressed at significantly higher levels with log (FC) > 2 and padj < 0.05. The RNA seq results were subjected to qRT-PCR validation. Among the nine genes selected, the expression of HAS1, TRIB2, P2RX1. PRG4, CNR2, and MMP25 were significantly higher (p < 0.05) in the anagen phase, whereas MC4R, GIPC2, and CDO1 were significantly expressed (p < 0.05) in the telogen phase which supports and validates the gene expression pattern from the RNA-sequencing. Differentially expressed genes revealed that Pashmina fiber initiation is largely controlled by signaling pathways like Wnt, NF-Kappa, JAK-STAT, Hippo, MAPK, Calcium, and PI3K-Akt. Expression of genes from the Integrin family, Cell adhesion molecules, and ECM-receptors were observed to be at much higher levels during anagen. We identified key genes (IL36RN, IGF2, ITGAV, ITGA5, ITCCR7, CXCL5, C3, CCL19, and CXCR3) and a collagen cluster which might be tightly correlated with anagen-induction. The regulatory network suggests the potential role of RUNX3, NR2F1/2, and GATA family transcription factors in anagen-initiation and maintaining fiber quality in Pashmina goats.

Pashmina goat (Capra hircus) is an economically important animal genetic resource adapted to very harsh cold arid agro-climatic conditions of Ladakh region (Jammu and Kashmir-India). The cold desert of Ladakh has a very short growing season and remains landlocked for more than half of the year. The mercury level of this landlocked high-altitude habitat (5500-6000 m above mean sea level) of Pashmina goat fluctuates between + 35 °C (short summers) and -40 °C (long winters). Under these stressful conditions (cold, arid, hypoxic and scanty vegetation), Pashmina goats remain active with different adaptation strategies. This goat produces the world's finest (11-14 µ) natural fibre (Pashmina fiber) which is used in the making of world-famous Pashmina/Cashmere shawls.
Pashmina fibre is the soft under-coat of the Pashmina goat mixed with coarse outer coat known as guard hair 1 . The guard hair develops from primary hair follicles (PHF) and Pashmina from secondary hair follicles (SHF). Pashmina fibre provides softness and luster to the product and has the glamour of being very rare. The unique www.nature.com/scientificreports/ fiber properties like fineness, texture and warmness are co-related with Pashmina genetics and breeding traits. The Pashmina fibers shed annually and are harvested by combing. Hair follicles (HF) undergo cycles of the rapid growth phase (anagen), regression phase (catagen) and nogrowth phase (telogen) throughput the lifetime. Hair follicles produce entire new hair shaft during the anagen phase. However, the underlying mechanism remains classical in the area 2 . Different approaches and animal models have been utilized to decipher the genetic basis of fiber transition [3][4][5][6][7][8][9] . However, asynchronous hair growth, the difference in skin anatomy and physiology restricts the identification of key molecular determinants mediating fiber transition.
Regulation of the hair cycle involves complex signaling interactions between Wnt (Wingless Integrin), Shh (Sonic HedgeHog), Notch, BMP (Bone Morphogenic Protein) and other signaling Pathways 6 . WNT 10,11 and Shh 12 signaling is indispensably important for new anagen, whereas BMPs 13 have been implicated in follicle differentiation. However, signaling pathways involved in hair follicle cycling are not sufficiently studied to date and new pathways are yet to be discovered for the proper understanding of the hair cycle 2 .
The secondary follicles of Pashmina goat present an excellent model for studying hair biology due to the circannual fiber cycle and synchronized fiber growth 14 . This follicle is also an excellent model for studying diverse cellular, molecular and biological processes 15,16 . In the present study transcriptomic profiling of skin biopsies containing secondary hair follicles were utilized to identify key genes and signaling pathways involved in hair follicle transition from no-growth phase (telogen) to growth phase (anagen).

Materials and methods
Experimental design. Changthangi goats were selected from the flock maintained at High Mountain Arid Agricultural Research Institute (HMAARI) Stakna Leh (Ladakh). All animals were kept under the identical conditions of natural photoperiod and natural temperatures. This study was approved by the Animal Welfare and Ethics Committee of the Sher-e-Kashmir University of Agricultural Science and Technology of Kashmir (SKUAST-K). All animal experiments were conducted in strict accordance with the rules and guidelines outlined by the SKUAST-K Animal Welfare Committee. As per the guidelines of the committee, the skin samples were collected aseptically using skin biopsy punch under local anaesthesia with minimal pain and discomfort to the animal.
Ten unrelated Pashmina goats of the same age (24 months) and sex (males) were repeatedly sampled during anagen (October active growth phase) and telogen (March-resting phase, before combing) 17 . The skin samples containing Pashmina follicles were collected from the flanking region of each goat. The skin biopsies samples were snap-frozen in liquid nitrogen and shipped to the laboratory in RNA-later for processing. All samples corresponding to a particular stage were collected at the same time.
Total RNA extraction, library construction, and sequencing. Total RNA was isolated from the skin tissues using the Trizol method (Invitrogen, USA) according to the manufacturer's protocol. RNA samples with a RIN value greater than 7.0 were selected for RNA-sequencing (RNA-seq). A total of 19 (10 anagen and 9 telogen) samples were further carried forward for the analysis as one of the samples from telogen did not qualify the minimum quality criteria for RNA-seq. For cDNA library preparation and sequencing, RNA samples were stored at − 80 °C. Approximately 4 µg of total RNA was used to prepare the RNA sequencing library using the TruSeq RNA Sample Prep Kits (Illumina) as per the kit's protocol. Agilent-tape station plots were used at every step to assess mRNA quality, enrichment success, fragmentation sizes, and final library sizes. Finally, the amplified fragments were sequenced using Illumina HiSeq 2500 to obtain 2 × 100 bp paired-end (PE) reads. Data analysis. The raw reads were pre-processed to remove the adapter sequences, low-quality reads and low-quality bases filtration towards 3′-end using cutadapt program v2.10 18 . Filtered reads were mapped to Capra hircus reference assembly ARS1 using Hisat program v2.20 (release date 2/6/2020) 19 . Quality control and alignment statistics (Supplementary Table 1) suggest the sequencing data were uniform among all sets of samples. Differential expression analysis between two contrast groups was performed using edgeR v3.30.3 20 using Trimmed Means of M-values (TMM) normalization method and paired test, after removing low expressed genes across all samples. DEGs between different cycling stages were screened based on the threshold of FDR corrected P-value < 0.05 and absolute log2 (fold change) > 1. Significantly, dysregulated genes were subjected to functional annotation and pathways enrichment analysis using KOBAS server v3.0 21,22 . Validation of DEGs with qPCR. cDNA was synthesized from 0.5 μg of the same total RNA used in RNAsequencing using the Revert Aid First Strand cDNA Synthesis Kit (Thermo Scientific, USA) as per the manufacturer's protocol. qPCR reactions were run on a Roche Lightcycler 480 II in a 20-μl reaction containing 0.5 μl of cDNA template, 10 μl of 2 × SYBR Green Master Mix, 0.3 μl of each primer (10 μmol/μl) and 8.9 μl nuclease-free water. The amplification program consisted of one cycle at 95 °C for 10 s, followed by 40 cycles of 95 °C for 15 s and 55 °C for 34 s. The qPCR reactions for each gene were performed with three biological replicates. Relative gene expression was normalized to the expression of goat GAPDH and calculated with the 2 −ΔΔCT method 23 . The expression levels of the genes obtained in RNA-seq and qPCR were compared with Pearson correlation coefficient.

Results
Primary processing of reads and differentially expressed genes. A  www.nature.com/scientificreports/ genomic regions. The characteristics of RNA-seq data were determined using PCA plots (Supplementary Fig. 1), showing a clear segregation and clustering of samples from different stages. Based on differential expression analysis using edgeR, 1094 genes were significantly dysregulated (log (FC) > 1, p adj < 0.05). Of the DEGs, a total of 748 were significantly up-regulated and 346 were down-regulated in the anagen (Supplementary Table 2) as compared to telogen. Two-dimensional hierarchical clustering and heatmap of top differentially expressed genes (q-value) were shown in Fig. 1.
The hair cycle is a highly regulated process, earlier studies suggest growth factors, cytokines, hormones, and transcription factors (TFs) play a critical role in mediating overall hair cycle 3 . In this study, a TF-regulatory network was generated (Fig. 3) to identify the critical TFs mediating anagen induction in Pashmina goats [25][26][27] . The |(log (FC))| and degree were used to identify key TFs in regulating hair cycle. The higher the two quantitative values of a gene, the more important it is in the TF-regulatory network. The TF-gene interaction network consisted of 50 nodes and 150 interactions (Fig. 3). The top-ranked TFs were Nuclear Receptor Subfamily 2 Group  Validation of DEGs with qPCR. The differential expression of 5 genes namely; HAS1, TRIB2, P2RX1, PRG4, CNR2, MMP25, GIPC2, CDO1, and MC4R was validated by qRT-PCR. Primer pairs for these genes are listed in Supplementary Table 4. The expression profile of these genes obtained by qRT-PCR showed a similar trend (Pearson's correlation coefficient = 0.88) with the RNA-seq results, thereby validating the RNA-seq results ( Supplementary Fig. 2).

Discussion and conclusion
Hair growth has been reported to be regulated by a complex mechanism involving multiple endogenous and exogenous factors. Understanding the genetic basis of fibre production phenotypes could contribute to the improvement of the fibre production efficiency in Pashmina goats and also help in identifying molecular determinants involved in anagen initiation. In this study, skin samples specifically containing SHFs were collected during anagen and telogen to identify key regulators and pathways associated with the telogen-anagen (TA) transition. Enrichment analysis of the DEGs for the TA transition revealed that the transition involves as many as 22 unique pathways. Among the enriched pathways Wnt, NF-B, JAK-STAT, MAPK and Calcium signaling are known in regulating hair follicle morphogenesis and development. Wnt acts as master regulating during hair follicle morphogenesis. NF-κB and Wnt signaling pathways play a vital role in hair follicle initiation and development. Strong NF-κB activity was detected in the secondary hair germ of late telogen and early anagen HFs in mouse 28 , suggesting a potential role for NF-κB in HF stem/progenitor cell activation during anagen induction. The present study also suggests a possible role of NF-κB in anagen induction. In the present study, the up-regulation of IL36RN during telogen which inhibits NF-κB suggests its potential role in maintaining telogen in Pashmina goats. The JAK-STAT pathway has recently been reported to be important in anagen induction and is responsible for jump-starting the hair cycle 29,30 . Up-regulation of IGF2 and enrichment of MAPK in anagen suggests the role of IGF2 in promoting anagen via MAPK signaling. Calcium signaling pathway was also significantly enriched. A transient role of Ca is reported during HF cycling 31 . Also, the melanogenesis pathway was significantly enriched in TA transition suggesting that melanogenesis pathway also resonate with hair cycle in Pashmina goats 32 .
The most distinguishing feature of the present study was significant enrichment of ECM-receptor interaction and Cell adhesion molecules (CAMs). Hair follicle growth depends on the interaction between CAMs and ECMreceptors (enriched by 18 genes and 24 genes respectively), it was found that ITGAV was strongly up-regulated in telogen samples, suggesting the potential role of ITGAV in maintaining hair follicle in the telogen phase. Protein-protein interaction network ( Fig. 2A) suggests the potential role of Integrin molecules like ITGB7, ITGA5, ITGA9 and ITGA11 in anagen-induction. In human adipose stem cells, the increase in ITGAV causes reduced cell proliferation 33 . However, the increased levels of ITGA5 activates cell proliferation and differentiation. Our data might postulate a functional role of ITGAV and ITGA5 in hair cycle.
A total of 10 key TFs were identified from the TF-regulatory network namely, GATA2, GATA1, RUNX3, NR2F1, NR2F2, HIC1, SPI1, MIXL1, ZFP2 and KLF1. Interestingly, ZFP2 was the only gene downregulated in the active growth phase (anagen). The GATA family of TFs plays an important role during embryonic development, including cell fate decision and tissue morphogenesis. GATA1 knockdown mouse exhibits developmental arrest and cell death 36,37 , network analysis suggest GATA1 interact with GATA2 may be for controlled cell differentiation and proliferation. RUNX3 plays a critical role in normal hair growth, a RUNX3 knockdown mice shows significant change in hair structure and composition 38 . This study also suggests the potential role of more diverse TFs (like NR2F1, NR2F2, HIC1, SPI1, MIXL1, ZFP2 and KLF1) in promoting anagen.

Data availability
The sequencing data is available in NCBI under accession number GSE164100.