Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Integrative analysis reveals ncRNA-mediated molecular regulatory network driving secondary hair follicle regression in cashmere goats

BMC Genomics volume 19, Article number: 222 (2018) Cite this article

Abstract

Background

Cashmere is a keratinized product derived from the secondary hair follicles (SHFs) of cashmere goat skins. The cashmere fiber stops growing following the transition from the actively proliferating anagen stage to the apoptosis-driven catagen stage. However, little is known regarding the molecular mechanisms responsible for the occurrence of apoptosis in SHFs, especially as pertains to the role of non-coding RNAs (ncRNAs) and their interactions with other molecules. Hair follicle (HF) degeneration is caused by localized apoptosis in the skin, while anti-apoptosis pathways may coexist in adjacent HFs. Thus, elucidating the molecular interactions responsible for apoptosis and anti-apoptosis in the skin will provide insights into HF regression.

Results

We used multiple-omics approaches to systematically identify long non-coding RNAs (lncRNAs), microRNAs (miRNAs) and mRNAs expressed in cashmere goat skins in two crucial phases (catagen vs. anagen) of HF growth. Skin samples were collected from three cashmere goats at the anagen (September) and catagen (February) stages, and six lncRNA libraries and six miRNA libraries were constructed for further analysis. We identified 1122 known and 403 novel lncRNAs in the goat skins, 173 of which were differentially expressed between the anagen and catagen stages. We further identified 3500 gene-encoding transcripts that were differentially expressed between these two phases. We also identified 411 known miRNAs and 307 novel miRNAs, including 72 differentially expressed miRNAs. We further investigated the target genes of lncRNAs via both cis- and trans-regulation during HF growth. Our data suggest that lncRNAs and miRNAs act synergistically in the HF growth transition, and the catagen inducer factors (TGFβ1 and BDNF) were regulated by miR-873 and lnc108635596 in the lncRNA-miRNA-mRNA networks.

Conclusion

This study enriches the repertoire of ncRNAs in goats and other mammals, and contributes to a better understanding of the molecular mechanisms of ncRNAs involved in the regulation of HF growth and regression in goats and other hair-producing species.

Background

Cashmere goats are the main livestock breed used for the production of both cashmere and meat [1]. Cashmere has an important status in the textile industry due to its high economic value [2]. Improvement the quantity and quality of cashmere becomes an important breeding goal, and molecular breeding is a convenient way to locate functional genes that increase fiber production [3, 4]. Preventing or delaying hair loss is a matter of some priority. Hair loss and production mechanisms remain incompletely explained, though there have been several studies involving hair regeneration in humans or other model animals [5, 6].

Cashmere is produced by the secondary hair follicles (SHFs), which exhibit an annual periodicity, undergoing anagen (growth), catagen (regression), and telogen (resting) phases annually. The majority of the time spent in each cycle is occupied by anagen, while catagen can have a critical effect on telogen or even the complete cycle. The basis for HF involution rests in the unique follicular epithelial and mesenchymal elements, as well as other cells type (adipocytes) intercellular molecules communication [7, 8]. Some of the molecular signals involved in HF regression process have been determined, including fibroblast growth factor (FGF), transforming growth factor-β (TGF-β), tumor necrosis factor-β (TNF-β), Wnt, sonic hedgehog (SHH), neurotrophins (NT), and homeobox proteins [7, 9,10,11,12,13].

HF regression is not only associated with the regulation of HF-structured cells but is also affected by other types of cellular changes in the skin environment. Whole-transcriptome sequencing could provide new insights into the molecular regulatory mechanisms of the HF cycle and the interactions among HF cells. The functions of microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) have been extensively reported in livestock, such as sheep and goat [14, 15], which played an important regulatory role in biological processes such as cell proliferation, differentiation, and apoptosis. MiRNAs are widely reported in animals and expressed in a temporal and spatial manner, playing an important role in hair follicle development and cycle [16,17,18,19]. We previously found that miRNAs are widely expressed in the skins of cashmere goats and that their expression levels were altered across the different HF cycling phases [16, 20, 21]. Specifically, the expression of miR-31 was significantly higher in the growth phase than in the regression phase; this miRNA targets the regulation of KRT16, KRT17, DLX3, and FGF10, thus affecting hair growth [22]. MiR-214 controls skin and HF development by modulating the activity of the Wnt pathway [23]. MiR-22 has been shown to be associated with HF degeneration and inhibits the expression of transcription factors such as DLX3, FOXN1, and HOXC13 [24]. Thus, miRNAs play an indispensable regulatory role in various biological processes during the HF cycle and in the HF transitions to other stages.

Other non-coding RNAs (ncRNAs), such as lncRNAs, are essential for the regulation of hair growth and the HF cycle, though the functions of the lncRNAs involved in the HF cycle remain unclear. The expression of lncRNAs in mouse dermal papilla cells (DPCs) changes with subsequent passage generations, indicating that lncRNAs are related to dermal papilla (DP) characteristics [25]. LncRNAs have been found to be associated with hair growth, playing an important role in the development of SHFs in sheep [15]. The lncRNA PlncRNA-1 regulates the proliferation and differentiation of HF stem cells through the TGFβ1-mediated Wnt/β-catenin signaling pathway [26]. The expression of lncRNA-H19 changes according to the growth phase of goat SHFs [27]. Overall, lncRNAs as well as miRNAs play an important role in the regulation of HF growth and development.

Despite this progress, the regulation of hair cycling in mammals is complex, and there may be other regulatory channels involved. Previous studies have reported that lncRNAs act as regulatory genes that compete with miRNAs [28] to not only directly inhibit mRNA expression but also bind miRNAs to regulate mRNA expression. In this study, we aimed to elucidate the molecular mechanism of HF regeneration by determining the expression levels of mRNAs, lncRNAs, and miRNAs and their corresponding relationships in the skin microenvironment.

Methods

Samples

Three two-year-old female Shanbei Cashmere goats with unrelated genetic background were used in this study. Skin samples were biopsied at mid-September and mid-February, as previously described [29]. To minimize animal suffering, procaine was used for local anesthesia. The goats were sampled from the Shanbei Cashmere Goat Farm of Hengshan, Yulin, China (located at 37°21′–38°14′ N and 108°56′–110°01′ E), being raised in the same environment. Dorsal skin samples were collected from between ribs 12 and 13. Each skin sample, about 1 cm2, was cut into pieces and then stored in an RNA/DNA sample protection reagent (Takara, Dalian, China), immediately. Samples were transported in dry ice and stored at − 80 °C for total RNA extraction. All sampling procedures in this experiment were in accordance with approved guidelines of the Animal Care and Use Committee of the Northwest A&F University (Approval ID: 2014ZX08008–002).

Total RNA isolation, library preparation, and sequencing

Total RNA was extracted using an Eastep® Super Total RNA Extraction Kit (Promega, Shanghai, China), according to the manufacturer’s instructions. We obtained two libraries from each sample: a lncRNA library and a miRNA library. The lncRNA library was prepared following a previous description [30], and library quality was assessed on the Agilent Bioanalyzer 2100 system. Libraries were sequenced on an Illumina HiSeq 4000 platform as 150-bp paired-end reads. Small RNA libraries were constructed and sequenced following a previous description [30], and the libraries were sequenced on an Illumina HiSeq 2500 platform using single-end reads.

Quality control, annotation, and expression levels

Raw reads from lncRNA libraries were routinely processed using a Perl script. In this step, adapters, reads containing over 10% Ns, and low-quality reads (> 50% of bases with Phred scores < 5%) were removed. The Phred score (Q20, Q30) and GC content of the clean data were calculated. All subsequent analyses were based on the high-quality data. The goat (Capra hircus) reference genome and gene annotation files (Ver. ASM170441v1) were downloaded from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/). An index of the goat reference genome was built using Bowtie2 v2.3.1 [31], and paired-end clean reads were aligned to the reference genome using TopHat2 v2.1.1 [32]. Cufflinks v2.1.1 [33] was used to analyze gene patterns. lncRNAs were identified based on their structures and the fact that they do not encode proteins; specifically, the following five criteria were used [34]: 1) transcripts must contain no fewer than 2 exons; 2) the transcript length must be more than 200 bp; 3) screening the known lncRNAs and the unknown transcript left for the following analysis; 4) quantification of each transcript must be no less than 0.5; 5) remaining unknown transcripts that potentially encoded proteins were assessed. Tophat2 was run with ‘–library-type fr-firststrand’, and Cufflinks was run with ‘min-frags-per-transfrag = 0’, while other parameters were set to the defaults.

Raw reads from small RNA libraries were first processed using custom Perl and Python scripts, and redundant regions were removed. Then, we chose a certain range of length from clean reads for all the downstream analyses, and the annotation was performed following previously described methods [35] with miRBase21.0 used as a reference. Cuffdiff (v2.1.1) [33] was also used to calculate the fragments per kb per million reads (FPKM) for both lncRNAs and coding genes in each sample. MiRNA expression levels were estimated by transcripts per million (TPM) as previously described [35]. For biological replication, transcripts, genes, or miRNAs with a P-value of < 0.05 were considered differentially expressed between the two groups of goats. DEseq2 was used to analyze all differential expression experiments.

Quantitative real-time PCR (qPCR) validation

Total RNAs were extracted from adult goat skin samples from the two groups and used for qPCR analysis. First strand cDNA was synthesized using the Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (#K1622, Thermo Scientific, USA) according to the manufacturer’s instructions and was then subjected to quantification using a standard SYBR Premix Ex Taq (Tli RNaseH Plus) kit (#RR420A, Takara, China) on the Bio-Rad CFX96 Real-Time PCR Detection System with β-actin as an endogenous control. The qPCR procedure was carried out according to the instructions of the reagent kit: pre-denaturation at 95 °C for 30 s, followed by 41 cycles of denaturation at 95 °C for 5 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s. The primers used for this experiment are listed in Additional file 1. Biological and technical replication was performed in triplicate for each sample. Relative gene expression was calculated using the 2-ΔΔCt method and quantified relative to β-actin. Student’s t-tests were used for statistical analysis, and a P-value < 0.05 was considered to be significant. Values are expressed as means ± SD, * P < 0.05, ** P < 0.01.

GO and KEGG pathway analysis

Gene ontology (GO) enrichment analysis of DE RNAs or lncRNA target genes was conducted with the GOseq R package, with a correction for gene length bias. GO terms with corrected P-values < 0.05 were considered significantly enriched DE genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) database is a resource for understanding the high-level functions and utilities of a biological system, such as a cell, organism, or ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg). We used KOBAS software to assess the statistical enrichment of DE genes or lncRNA target genes in KEGG pathways. Directed acyclic graph (DAG) is a graphical display of DE genes GO enrichment analysis results, and DAGs were drawn in the biological process (BP), cellular component (CC) and molecular function (MF), respectively.

Correlation and co-expression/co-localization analysis

Co-expression analysis was based on calculating the Pearson correlation coefficients (PCCs) between coding genes and noncoding transcripts according to their expression levels. An absolute value of the parameter PCC ≥ 0.95, P-value < 0.01, and false discovery rate (FDR) < 0.01 was used for identifying genes for further analysis. For identifying cis-regulation, lncRNAs that act on neighboring target genes were investigated. We searched for coding genes 10 k/100 k upstream and downstream of each lncRNA and then analyzed their functions. For identifying trans-regulation, lncRNAs and target genes were identified based on their expression levels. We calculated the expression level correlation between lncRNAs and coding genes using custom scripts.

Competing endogenous RNAs (CeRNAs) network analysis

To reveal the roles and interactions of ncRNAs and mRNAs in the HF growth cycle, we constructed ncRNA regulatory networks. Those lncRNAs and mRNAs whose expression levels were meaningfully correlated were included in this analysis. Potential miRNA response elements were searched for the sequences of lncRNAs and mRNAs, and we identified overlaps in predicted miRNA seed sequence binding sites and lncRNAs binding sites in the target mRNA as part of the lncRNA–miRNA–mRNA interaction. The miRNA binding sites were predicted by miRBase (http://mirbase.orghttp://mirbase.org), while the miRNA–mRNA interactions were predicted by Targetscan (http://www.targetscan.org/). The interaction network was built and visually displayed using Cytoscape software based on the screening of lncRNA–miRNA–mRNA pairs.

Results

Morphological features of goat skins at different stages

To examine the morphological differences in the growth to regression stages, we collected cashmere goat skin samples at anagen (September) and catagen (February). These time points were determined according to HF morphological analysis and fiber growth features. Hematoxylin and eosin (H&E) staining indicated clear differences in the hair bulbs of anagen and catagen goat skins. Specially, the sharp/size of DPs had a great change and their appearance were atrophying in February, while they were plump in September (Fig. 1a). These results indicate that a series of biological processes occur in goat skins throughout the year, and this further guaranteed that the proper samples were collected for further analysis. Therefore, a comprehensive whole-transcriptome sequencing analysis was used to systematically investigate the role of genes and ncRNAs in the biological processes of the HF transition from the anagen to the catagen stage (Fig. 1b).

Fig. 1
figure 1

H&E staining of goat skins and overall study design. (a) H&E staining results of goat skins in the anagen and catagen phases. Red arrows show the location of hair dermal papilla. (b) Schematic workflow of the experimental design of this study

Full size image

Identification of lncRNAs and miRNAs in goat skin

An average of 183,914,710 raw reads were produced on the Illumina HiSeq 4000 platform for each sample. After filtering out adaptor-related and low-quality sequences and those containing Ns, we obtained 176,158,389 clean reads, accounting for 26.4 Gb of clean sequence for each sample. Subsequently, we classified and mapped the clean reads to the latest goat reference genome assembly (ARS1, https://www.ncbi.nlm.nih.gov/assembly/GCF_001704415.1). The output data, raw read classification, and mapping region for each sample are summarized in Additional file 2 (Tables S1–S2) and Additional file 3: Figure S1, Additional file 3: Figure S2). We identified lncRNAs by characteristics that distinguish them from other RNAs (coding mRNAs, tRNAs, rRNAs, snRNAs, snoRNAs, pre-miRNAs, and pseudogenes), including the use of five rigorous criteria and screening tools for identifying potential coding sequences (see Additional file 3: Figure S3a). We identified 45,638 transcripts, 1122 known lncRNAs, 403 novel lncRNAs (Additional file 3: Figure S3b), and 350 transcripts of uncertain coding potential (TUCP). A list of these and their expression levels is provided in Additional file 4. For miRNA-seq, we obtained more than 28.6 M clean reads, accounting for 1.43 Gb of clean sequence for each sample (Additional file 2: Table S3). We identified 411 known miRNAs and 307 novel miRNAs (Additional file 5). A summary of the identified RNAs is provided in Table 1, and those determined to be differentially expressed (DE) between the anagen and catagen stages were used for further analysis.

Table 1 Summary of identified genes and ncRNAs
Full size table

Differentially expressed mRNA, lncRNA, and miRNA profiles

To determine whether ncRNAs are involved in the HF regression process. DE ncRNAs and mRNAs from catagen and anagen stages were visualized using a volcano plot and clustering map, while overlapping RNAs expressed in the two groups were visualized using a Venn diagram. There were 173 DE lncRNAs comparing the catagen to the anagen, including 98 up-regulated and 75 down-regulated lncRNAs (Fig. 2a–b, Additional file 6). Similarly, there were 3500 DE mRNA transcripts (3357 genes) comparing the catagen to the anagen, including 1830 up-regulated and 1670 down-regulated transcripts (Fig. 3a–c, Additional file 7), and there were 72 DE miRNAs comparing the catagen to the anagen, including 39 up-regulated and 33 down-regulated miRNAs (Fig. 4a–c, Additional file 8). Fig. 5a shows a summary histogram of DE lncRNAs, mRNAs, and miRNAs, and the top 20 most significantly DE ncRNAs are provided in Table 2. We found that those lncRNAs (LOC102190274, LOC108635596, LOC108635657, LOC108636746, LOC108635658, LOC102188339, LOC108635659, and LOC108635656) are potential regulators. LncRNAs regulate the expression of target genes (mRNAs), and this can be demonstrated by co-localization and co-expression. Given that DE mRNAs could be directly or indirectly regulated by lncRNAs, we identified the overlap between lncRNA target genes and DE mRNAs for further analysis. Fig. 5b shows the target mRNAs of lncRNAs, based on co-localization and co-expression, and the overlap between these and DE mRNAs in a Venn diagram.

Fig. 2
figure 2

LncRNA expression profile changes in goat skins. (a) Volcano plot indicating up- and down-regulated lncRNAs in the catagen stage compared with the anagen stage. (b) Heat map of lncRNAs showing hierarchical clustering of altered lncRNAs in the catagen stage compared with the anagen stage. Up- and down-regulated genes are in red and green, respectively. (c) Venn diagram showing the number of overlapping lncRNAs in the catagen and anagen stages

Full size image
Fig. 3
figure 3

mRNA expression profile changes in goat skins. (a) Volcano plot indicating up- and down-regulated mRNA transcripts in the catagen stage compared with the anagen stage. (b) Heat map of mRNA transcripts showing hierarchical clustering of altered mRNA transcripts in the catagen stage compared with the anagen stage. Up- and down-regulated genes are in red and green, respectively. (c) Venn diagram showing the number of overlapping mRNA transcripts in the catagen and anagen stages

Full size image
Fig. 4
figure 4

MiRNA expression profile changes in goat skins. (a) Volcano plot indicating up- and down-regulated miRNAs in the catagen stage compared with the anagen stage. (b) Heat map of mRNA transcripts showing hierarchical clustering of altered miRNAs in the catagen stage compared with the anagen stage. Up- and down-regulated genes are in red and green, respectively. (c) Venn diagram showing the number of overlapping miRNAs in the catagen and anagen stages

Full size image
Fig. 5
figure 5

Count of relative ncRNAs and mRNAs in goat skins. (a) Histogram showing the number of up- and down-regulated ncRNAs and miRNAs in goat skins. (b) Venn diagram showing the number of overlapping targeted mRNAs of up-regulated lncRNAs, targeted mRNAs of down-regulated lncRNAs, and up- and down-regulated mRNAs

Full size image
Table 2 The information of top 20 ncRNAs
Full size table

In comparing the RNAs expressed in the catagen and anagen phases, there were some genes of note. As keratin protein (KRT) and keratin-associated protein (KRTAP) are strongly associated with hair growth, we focused on these. We found that 74 KRT and KRTAP transcripts were differentially expressed between the catagen and the anagen stages (Fig. 6a). These genes are involved in driving the physiological characteristics of hair growth. KRT and KRTAP comprise the hair shaft and are an important part of cashmere. Most DE KRT and KRTAP genes (KRT38, KRT4, KRTAP15–1, KRTAP13.1, and KRTAP3–1) are involved in the construction of hair and were more highly expressed in the anagen than in the catagen stage. We also found that a few genes (KRT2, KRTDAP, KRT77, and KRT80) involved in epidermis keratinization were more highly expressed in the catagen stage. We also investigated DE growth factors or ligands associated with hair growth, and we found that the expression of growth arrest-specific genes (GAS1, GAS6, and GAS7) was higher in the catagen than in the anagen stage. We also found that growth factors and their receptors (FGF22, FGF21, FGF2, GDF11, IGF1, and FGFR4), as well as up-regulated skeletal muscle growth 5 homolog (USMG5) were more highly expressed in the anagen stage (Fig. 6b).

Fig. 6
figure 6

Heat maps for gene groups. (a) Heat map of differentially expressed KRT and KRTAP genes. (b) Heat map of differentially expressed growth factor-related genes

Full size image

Validation of mRNAs and lncRNAs

To validate the reliability of the sequencing results, the expression changes of some mRNAs (WNT11, PLIN4, LAMA5, LAMA2, PLIN1, MSX2, HSP70.1, TGFBR2, and WNT4) and lncRNAs (lnc10218839, lnc102190274, lnc102189880, lnc106502925, and lnc102171315) were validated by qPCR (Fig. 7a). The expression levels of the genes as determined by the sequencing and qPCR methods are shown in Fig. 7b. These exhibited a correlation coefficient of 0.865 and P-value of 6.599E-05. The qPCR expression levels of all validated mRNAs and lncRNAs were consistent with the results obtained from HiSeq.

Fig. 7
figure 7

qPCR validation of RNA-seq. (a) Nine coding genes and five noncoding genes (lncRNAs are renamed “lnc” + the gene accession number). (b) Correlation analysis between qPCR and RNA-seq results. The x-axis indicates the fold change value log2ratio(catagen/anagen) according to qPCR, and the y-axis indicates the log2ratio(FPKMcatagen/anagen) value according to RNA-seq

Full size image

Functional prediction of ncRNAs in the HF growth cycle

To explore the potential regulatory roles of lncRNAs in the HF growth cycle, we performed an integrated co-expression and co-localization network analysis of DE lncRNAs and mRNAs. We further selected genes in which the absolute value of the correlation was > 0.95 to predict the functions of lncRNAs using GO and KEGG analysis tools. Additional file 3: Figure S4a–c (Additional file 3) shows the BP, CC, and MF categories of the GO enrichment analysis based on the co-expression and co-localization of DE lncRNAs. The most significantly enriched BP terms were cellular metabolic process, metabolic process, and intracellular transport, and the noteworthy CC terms were protein complex, intracellular organelle part, and intermediate filament. The most significantly enriched MF terms were binding, protein binding, and ion binding (Additional file 3: Figure S4d). Additional file 3:Figure S5a–c (Additional file 3) shows the BP, CC, and MF categories of the GO enrichment analysis involving the overlapping co-expressed and co-localized DE lncRNA target genes, DE miRNA target genes, and predicted mRNAs. The most significantly enriched BP terms were biosynthetic process, metabolic process, and gene expression, and the most significantly enriched CC terms were intermediate filament, cytoskeleton, and protein complex. The most significantly enriched MF terms were binding, catalytic activity, and cytokine activity (Additional file 3: Figure S5d). Additional file 3: Figure S6a–c (Additional file 3) shown the BP, CC, and MF categories of the GO enrichment analysis of all DE genes. The most significantly enriched BP terms were protein import into nucleus, protein localization to nucleus, and protein targeting to nucleus, while the most significantly enriched CC terms were intermediate filament, intermediate filament cytoskeleton, and keratin filament. The most significantly enriched MF terms were cytokine activity, chemokine activity, and binding (Additional file 3: Figure S6d). These enriched categories of DE ncRNAs and their target genes showed potential values and will be the focus of future studies in HF regression.

An enriched scatter diagram of the candidate genes provides a graphic display of the KEGG enrichment analysis. The degree of KEGG enrichment is assessed by the richness factor, Q-value, and number of genes. We listed the top 20 most enriched pathways for DE mRNAs and lncRNA-miRNA-mRNAs by Q-value range (Additional file 3: Figure S7a–b). The most enriched pathways were disease pathways (pathways in cancer and systemic lupus erythematosus) and signaling molecules and interaction pathways (extracellular matrix (ECM)-receptor interaction). The ncRNA target genes exhibited KEGG enrichment patterns similar to those of the DE mRNAs. We divided the DE mRNA KEGG enrichment pathways into two categories: up-regulated and down-regulated (Additional file 3: Figure S7c–d). In the up-regulated category, the most significant pathways were pathways in cancer, cytokine–cytokine receptor interactions, adherens junctions, and the TNF signaling pathway. Most of these participate in cell apoptosis. Among the down-regulated pathways, the most significant pathways were metabolic pathways, oxidative phosphorylation, cell cycle, and ECM-receptor interactions. Additional KEGG pathways are provided in Additional file 9, and they deserve further analysis in the future.

lncRNA-miRNA-mRNA networks

To explore the molecular mechanisms by which ncRNAs are involved in HF development and the HF cycle, we performed regulatory network analysis of ncRNAs and mRNAs in the two growth phases. MiRNAs have broadly regulatory functions involving RNA silencing and post-transcriptional regulation of gene expression, while lncRNAs also have extensive regulatory functions. We performed an integrated ncRNA and mRNA profiling analysis on the basis of their binding sites, with lncRNA-gene pairs and miRNA-gene pairs that share the same binding site being considered ceRNAs. We constructed lncRNA-miRNA-mRNA groups with the lncRNA as the decoy, miRNA as the center, and mRNA as the target (Fig. 8a). We then divided the ceRNA regulatory networks into two groups based on their gene expression patterns: up-down-up (Fig. 8b) and down-up-down (Fig. 8c). The results indicated that the expression of genes in the skin during HF regression is regulated by an ncRNA regulatory network that involves the action of ceRNAs. Our results (Table 3) present the regulatory relationships between ncRNAs and mRNAs in the processes of HF development and the HF cycle. Thus, further studies should be undertaken to better understand the mechanisms of HF development and cycle in the future.

Fig. 8
figure 8

LncRNA–miRNA–mRNA regulatory networks in goat skin. Circles indicate miRNAs, squares indicate coding genes, and triangles indicate lncRNAs. Red indicates up-regulation, and green indicates down-regulation (catagen/anagen)

Full size image
Table 3 ncRNAs and their potential target genes involved in HF cycle
Full size table

Discussion

Skin, as the first line of defense consisting of diverse cells (adipocyte and fibroblast) and affiliate organs (HFs and sweat glands), covers the entire body. HFs are an important component of the skin, consisting of primary HFs and SHFs. SHFs account for over 90% of the two types of HFs in cashmere goats, and they produce fine fibers in a rigorous annual cycle that involves anagen, catagen, and telogen [36]. During catagen, the DPs of SHFs shrink, and therefore, it is presumed that there is strong apoptosis signaling in DPs and anti-apoptosis singling in adjacent DPs. While seven theories have been proposed to explain the HF cycle (epithelial theory, papilla morphogen theory, bulge activation theory, resonance theory, oscillating signal theory, inherent embryonic cycle theory, and inhibition-disinhibition theory) [37], the underlying mechanism for regression remains unknown. In this study, we performed an ncRNA and coding RNA profiling analysis by comparing the DE RNAs of the anagen and the catagen phases, providing novel insights into the mechanism of the HF cycle.

We initially performed a histological analysis to show that the DPs of SHFs experienced atrophy in February and were full in September. This phenomenon was consistent with expectations and provided the basis for the subsequent analysis. Next, we performed whole-transcriptome sequencing and a comparative analysis of goat skins in different growth phases. Phenotypic appearance of HFs is affected by the expression of genes, which is regulated by various factors [38, 39]. The findings obtained from several studies have revealed the critical roles of ncRNAs in the hair growth cycle. Indeed, miRNAs have been frequently reported as regulators of skin and HF development [40], and there is also considerable evidence supporting the idea that lncRNAs play causal roles in the HF growth cycle [15, 25]. However, there was previously a lack of information regarding the role of ncRNAs in the HF cycle. In this study, we found that abundant RNAs, including lncRNAs, miRNAs, and coding mRNAs, were enriched in the skin and differentially expressed between the anagen and catagen phases. Some of these transcripts were differentially expressed in the skin at different growth phases, indicating that they play an important role in the regulation of biological processes, in particular those associated with apoptosis, anti-apoptosis, and growth processes [41,42,43]. Our results showed that a total of 173 lncRNAs, 72 miRNAs, and 3500 mRNA transcripts were significantly differentially regulated when comparing the catagen and anagen phases. These data indicate that many miRNAs are involved in HF development, and some of these miRNAs differ from those observed in previous research studies [16], although there were also some common miRNAs (miR-27a-5p, miR-99a-5p, miR-380-3p, and miR-9-5p) that deserve further attention. One reason for this discrepancy is the difference in analytical strategies. Another reason is that the same biological processes may be induced by a different combination of genes that coordinate with each other. In addition, the analysis of lncRNAs and their roles as ceRNAs provide new insight into these differences. Despite the fact that the functions of lncRNAs in HF development and the HF cycle are still poorly understood, the rhythmic expression of lncRNA genes indicates their casual roles in biological processes. Numerous lncRNAs, acting as decoys for miRNAs, have been found to play critical roles in biological functions such as apoptosis, proliferation, and differentiation [44,45,46,47].

Although we screened for DE ncRNAs from goat skins in different growth phases and a few ncRNAs were confirmed in the present study, the underlying mechanisms of ncRNAs in the HF cycle are poorly understood. We then analyzed DE ncRNAs-related gene functions and their corresponding pathways through GO and KEGG term enrichment analyses. Our data showed that the most significantly enriched pathway was pathways in cancer, a comprehensive pathway involving multiple cellular processes. Other cancer pathways, such as prostate cancer and proteoglycans in cancer, as well as cytokine–cytokine receptor interaction, estrogen signaling pathway, adherens junction, and the Jak-STAT signaling pathway pointed to roles in cell growth, proliferation, and apoptosis/anti-apoptosis. These pathways fully illustrate the balance between apoptosis and anti-apoptosis in the skin. The DPs of SHFs begin to atrophy while other cells experience accelerated growth and production of keratins to protect against environmental stress. There were also some lipid metabolism and endocrine system KEGG pathways contributing to fat deposition as well as the hair growth cycle. The predicted functions of ncRNAs in the HF cycle, as determined by the GO and KEGG analyses, should be studied at a greater depth in future work.

The majority of fine hair growth slows down or ceases during catagen, including the production of KRT and KRTAP, which are indispensable structural components of hair. Most of DE KRT and KRTAP genes were down-regulated during catagen. Interestingly, we found that the expression levels of KRT2, KRTDAP, KRT77, and KRT80 were higher during catagen. These genes mainly participate in epidermis keratinization, not in the production of hair components, indicating that epidermal cells experience extra growth during the catagen phase; this may be a defense mechanism against upcoming extreme environments (i.e., the cold winter). We also found that the expression of GAS1 was higher during catagen than during anagen; this gene plays a role in growth suppression, blocking entry to S phase and preventing the cycling of cells [48]. This is consistent with what happens to the DPs of SHFs during catagen. Additionally, the expression levels of FGF22, FGF21, FGF2, GDF11, IGF1, FGFR4, and USMG5 were lower during catagen, indicating that the growth of the entire skin slowed down, perhaps caused by the transition of SHFs from anagen to catagen. It is worth mentioning that HSPH1, HSPA6, HSP70.1, LOC102178315 (heat shock 70 kDa protein 1B), and HSPB1 were also highly expressed during catagen. These genes belong to the heat shock protein family and are involved in stress resistance, such as ensuring the correct folding of proteins or controlling the targeting of proteins for subsequent degradation [49,50,51]. These proteins would therefore be useful in preventing apoptosis, making a response to DPs adjacent organization. Next, we found that the expression levels of AR and DKK1 were higher during catagen. AR, an androgen receptor, modifies the expression of the Wnt antagonist DKK1 in DPs, preventing HF stem cell differentiation [52]. Thus, this is another factor associated with fatty acids or adipocytokines that might trigger a non-apoptotic DP cell death. Intradermal adipocytes are distributed around SHFs in the skin, making this hypothesis more relevant. We detected that the expression levels of FASN and TRL4 were higher during catagen, and these two genes may also play a role in non-apoptotic cell death [53].

Based on these lines of evidence, we propose a model (Fig. 9) for understanding what occurs in the skin during catagen, with gene locations based on previous studies [54]. Some evidence indicates that BMP2, TGFβ1, and BDNF could represent inducers of catagen, while IGF1 could inhibit catagen [23]. HF stem cells return to quiescence via the action of BMPs (BMP2) and Wnt inhibitors, and their characteristics are maintained during HF regression by LHX2, SOX2, ITGA3, ITGB4, and other factors [55]. Meanwhile, the matrix undergoes apoptosis, and hair shaft formation slows. HOXC13, MSX1, MSX2, and other factors regulate the proliferation and differentiation of the hair matrix, and their expression were down-regulated during anagen in our study. Under the mediation of various genes, apoptosis and anti-apoptosis coexist among different cells in a single organism. The expression levels of CASP8 and MYC are higher during catagen, contributing to apoptosis, while the increased expression levels of BAG3, HSP70, and XIAP during catagen may contribute to anti-apoptosis. CASP14 is helpful for keratinization. The factors leading to the regression of HF still require further study and additional lines of evidence.

Fig. 9
figure 9

A proposed model for the microenvironment of the skin during the HF transition from anagen to catagen. Red squares indicate that expression of these genes is up-regulated during catagen, while green squares indicate down-regulation

Full size image

Conclusions

In this study, we provided new insight into the molecular mechanisms driving the transition from anagen to catagen in HFs. We integrated analysis of lncRNAs, miRNAs, and mRNAs to determine their predicted roles and regulatory relationships in the HF cycle. We also provided a catalog of predicted ncRNAs in goat skin that will help to explain their regulatory roles in the goat HF cycle. In future studies, we plan to further examine the functions of these predicted DE ncRNAs, which may ultimately enable the full elucidation of the regulatory mechanisms associated with the goat HF cycle at the molecular level.

Abbreviations

BP:

Biological process

CC:

Cellular component

ceRNAs:

Competing endogenous RNAs

DEGs:

Differentially expressed genes

DP:

Dermal papilla

DPCs:

Dermal papilla cells

FPKM:

Fragments Per Kilobase of exon per Million fragments mapped

GO:

Gene Ontology

HF:

Hair follicle

KEGG:

Kyoto Encyclopedia of Genes and Genomes.

KRT:

Keratin protein

KRTAP:

Keratin associated protein

LncRNA:

Long noncoding RNA; miRNA, microRNA

MF:

Molecular function

SHF:

Secondary hair follicle

TMP:

Transcript per million

TUCP:

Transcripts of uncertain coding potential

USMG5:

Skeletal muscle growth 5 homolog

References

  1. Ansari-Renani HR. Cashmere production, harvesting, marketing and processing by nomads of Iran. Pastoralism. 2015;5(1):18.

    Article  Google Scholar 

  2. Waldron S, Brown C, Komarek AM. The Chinese cashmere industry: a global value chain analysis. Dev Policy Rev. 2014;32(5):589–610.

    Article  Google Scholar 

  3. Jin M, Wang L, Li S, Xing MX, Zhang X. Characterization and expression analysis of KAP7.1, KAP8.2 gene in Liaoning new-breeding cashmere goat hair follicle. Mol Biol Rep. 2011;38(5):3023–8.

    Article  CAS  PubMed  Google Scholar 

  4. Xu T, Guo XD, Wang H, Du XY, Gao XY, Liu DJ. De novo transcriptome assembly and differential gene expression profiling of three Capra hircus skin types during Anagen of the hair growth cycle. Int J Genomics. 2013;2013(10):269191.

    PubMed  PubMed Central  Google Scholar 

  5. Ramos PM, Miot HA. Female pattern hair loss: a clinical and pathophysiological review. An Bras Dermatol. 2015;90(4):529–43.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Plikus MV, Mayer JA, de la Cruz D, Baker RE, Maini PK, Maxson R, Chuong CM. Cyclic dermal BMP signalling regulates stem cell activation during hair regeneration. Nature. 2008;451(7176):340–U348.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Sennett R, Rendl M. Mesenchymal-epithelial interactions during hair follicle morphogenesis and cycling. Semin Cell Dev Biol. 2012;23(8):917–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Schmidt B, Horsley V. Unravelling hair follicle-adipocyte communication. Exp Dermatol. 2012;21(11):827–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Lee J, Tumbar T. Hairy tale of signaling in hair follicle development and cycling. Semin Cell Dev Biol. 2012;23(8):906–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Tsai SY, Sennett R, Rezza A, Clavel C, Grisanti L, Zemla R, Najam S, Rendl M. Wnt/beta-catenin signaling in dermal condensates is required for hair follicle formation. Dev Biol. 2014;385(2):179–88.

    Article  CAS  PubMed  Google Scholar 

  11. Qiu WM, Lei MX, Tang H, Yan HT, Wen XH, Zhang W, Tan RJ, Wang D, Wu JJ. Hoxc13 is a crucial regulator of murine hair cycle. Cell Tissue Res. 2016;364(1):149–58.

    Article  CAS  PubMed  Google Scholar 

  12. Rishikaysh P, Dev K, Diaz D, Qureshi WMS, Filip S, Mokry J. Signaling involved in hair follicle morphogenesis and development. Int J Mol Sci. 2014;15(1):1647–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Lindner G, Botchkarev VA, Botchkareva NV, Ling G, vanderVeen C, Paus R. Analysis of apoptosis during hair follicle regression (catagen). Am J Pathol. 1997;151(6):1601–17.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang S, Ge W, Luo Z, Guo Y, Jiao B, Qu L, Zhang Z, Wang X. Integrated analysis of coding genes and non-coding RNAs during hair follicle cycle of cashmere goat (Capra hircus). BMC Genomics. 2017;18(1):767.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Yue Y, Guo T, Yuan C, Liu J, Guo J, Feng R, Niu C, Sun X, Yang B. Integrated analysis of the roles of long noncoding RNA and coding RNA expression in sheep (Ovis aries) skin during initiation of secondary hair follicle. PLoS One. 2016;11(6):e0156890.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Yuan C, Wang X, Geng R, He X, Qu L, Chen Y. Discovery of cashmere goat (Capra hircus) microRNAs in skin and hair follicles by Solexa sequencing. BMC Genomics. 2013;14:511.

  17. Riccardi S, Bergling S, Sigoillot F, Beibel M, Werner A, Leighton-Davies J, Knehr J, Bouwmeester T, Parker CN, Roma G, et al. MiR-210 promotes sensory hair cell formation in the organ of corti. BMC Genomics. 2016;17:309.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Liu G, Liu R, Li Q, Tang X, Yu M, Li X, Cao J, Zhao S. Identification of microRNAs in wool follicles during anagen, catagen, and telogen phases in Tibetan sheep. PLoS One. 2013;8(10):e77801.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Chen X, Ge K, Wang M, Zhang C, Geng Z. Integrative analysis of the Pekin duck (Anas anas) MicroRNAome during feather follicle development. BMC Dev Biol. 2017;17(1):12.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Liu ZH, Xiao HM, Li HP, Zhao YH, Lai SY, Yu XL, Cai T, Du CG, Zhang WG, Li JQ. Identification of conserved and novel microRNAs in cashmere goat skin by deep sequencing. PLoS One. 2012;7(12):e50001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Wu ZY, Fu YH, Cao JH, Yu M, Tang XH, Zhao SH. Identification of differentially expressed miRNAs between white and black hair follicles by RNA-sequencing in the goat (Capra hircus). Int J Mol Sci. 2014;15(6):9531–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Mardaryev AN, Ahmed MI, Vlahov NV, Fessing MY, Gill JH, Sharov AA, Botchkareva NV. Micro-RNA-31 controls hair cycle-associated changes in gene expression programs of the skin and hair follicle. FASEB J. 2010;24(10):3869–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Ahmed MI, Alam M, Emelianov VU, Poterlowicz K, Patel A, Sharov AA, Mardaryev AN, Botchkareva NV. MicroRNA-214 controls skin and hair follicle development by modulating the activity of the Wnt pathway. J Cell Biol. 2014;207(4):549–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Yuan SK, Li FF, Meng QY, Zhao YQ, Chen L, Zhang HQ, Xue LX, Zhang XQ, Lengner C, Yu ZQ. Post-transcriptional regulation of keratinocyte progenitor cell expansion, differentiation and hair follicle regression by miR-22. PLoS Genet. 2015;11(5):e1005253.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Lin CM, Liu Y, Huang K, Chen XC, Cai BZ, Li HH, Yuan YP, Zhang H, Li Y. Long noncoding RNA expression in dermal papilla cells contributes to hairy gene regulation. Biochem Biophys Res Commun. 2014;453(3):508–14.

    Article  CAS  PubMed  Google Scholar 

  26. Si Y, Bai J, Wu J, Li Q, Mo Y, Fang R, Lai W. LncRNA PlncRNA1 regulates proliferation and differentiation of hair follicle stem cells through TGFbeta1mediated Wnt/betacatenin signal pathway. Mol Med Rep. 2017;17(1):1191–7.

    PubMed  Google Scholar 

  27. Zhu YB, Wang ZY, Yin RH, Jiao Q, Zhao SJ, Cong YY, Xue HL, Guo D, Wang SQ, Zhu YX, et al. A lncRNA-H19 transcript from secondary hair follicle of Liaoning cashmere goat: identification, regulatory network and expression regulated potentially by its promoter methylation. Gene. 2018;641:78–85.

    Article  CAS  PubMed  Google Scholar 

  28. Yoon JH, Abdelmohsen K, Gorospe M. Posttranscriptional gene regulation by long noncoding RNA. J Mol Biol. 2013;425(19):3723–30.

    Article  CAS  PubMed  Google Scholar 

  29. Gao Y, Wang XL, Yan HL, Zeng J, Ma S, Niu YY, Zhou GX, Jiang Y, Chen YL. Comparative transcriptome analysis of fetal skin reveals key genes related to hair follicle morphogenesis in cashmere goats. PLoS One. 2016;11(3):e0151118.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Ren H, Wang G, Chen L, Jiang J, Liu L, Li N, Zhao J, Sun X, Zhou P. Genome-wide analysis of long non-coding RNAs at early stage of skin pigmentation in goats (Capra hircus). BMC Genomics. 2016;17:67.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011;25(18):1915–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, Zhao X, Liang C, Wang Y, Sun L, et al. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One. 2010;5(12):e15224.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Shi B, Ding Q, He X, Zhu H, Niu Y, Cai B, Cai J, Lei A, Kang D, Yan H, et al. Tbeta4-overexpression based on the piggyBac transposon system in cashmere goats alters hair fiber characteristics. Transgenic Res. 2017;26(1):77–85.

    Article  CAS  PubMed  Google Scholar 

  37. Stenn KS, Paus R. Controls of hair follicle cycling. Physiol Rev. 2001;81(1):449–94.

    Article  CAS  PubMed  Google Scholar 

  38. Chi W, Wu E, Morgan BA. Earlier-born secondary hair follicles exhibit phenotypic plasticity. Exp Dermatol. 2015;24(4):265–8.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Yano K, Brown LF, Detmar M. Control of hair growth and follicle size by VEGF-mediated angiogenesis. J Clin Invest. 2001;107(4):409–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Cai T, Liu ZH, Wang ZX, Zhao M, Ju HL, Li JQ. miRNA in regulation of skin and hair follicle development. Yi Chuan. 2013;35(9):1087–94.

    Article  CAS  PubMed  Google Scholar 

  41. Pileczki V, Cojocneanu-Petric R, Maralani M, Neagoe IB, Sandulescu R. MicroRNAs as regulators of apoptosis mechanisms in cancer. Clujul Med. 2016;89(1):50–5.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Zhang Y, Takahashi S, Tasaka A, Yoshima T, Ochi H, Chayama K. Involvement of microRNA-224 in cell proliferation, migration, invasion, and anti-apoptosis in hepatocellular carcinoma. J Gastroenterol Hepatol. 2013;28(3):565–75.

    Article  CAS  PubMed  Google Scholar 

  43. Hwang HW, Mendell JT. MicroRNAs in cell proliferation, cell death, and tumorigenesis. Br J Cancer. 2007;96(Suppl):R40–4.

    PubMed  Google Scholar 

  44. Chen SC, Wang MD, Yang H, Mao L, He QW, Jin HJ, Ye ZM, Luo XY, Xia YP, Hu B. LncRNA TUG1 sponges microRNA-9 to promote neurons apoptosis by up-regulated Bcl2l11 under ischemia. Biochem Bioph Res Co. 2017;485(1):167–73.

    Article  Google Scholar 

  45. Hu K, Zhang J, Liang M. LncRNA AK015322 promotes proliferation of spermatogonial stem cell C18-4 by acting as a decoy for microRNA-19b-3p. In Vitro Cell Dev Biol Anim. 2017;53(3):277–84.

    Article  CAS  PubMed  Google Scholar 

  46. Pickard MR, Mourtada-Maarabouni M, Williams GT. Long non-coding RNA GAS5 regulates apoptosis in prostate cancer cell lines. Bba-Mol Basis Dis. 2013;1832(10):1613–23.

    Article  CAS  Google Scholar 

  47. Yu C, Li L, Xie F, Guo S, Liu F, Dong N, Wang Y. LncRNA TUG1 sponges miR-204-5p to promote osteoblast differentiation through upregulating Runx2 in aortic valve calcification. Cardiovasc Res. 2017;114(1):168–79.

    Article  Google Scholar 

  48. Li QG, Qin Y, Wei P, Lian P, Li YQ, Xu Y, Li XX, Li DW, Cai SJ. Gas1 inhibits metastatic and metabolic phenotypes in colorectal carcinoma. Mol Cancer Res. 2016;14(9):830–40.

    Article  CAS  PubMed  Google Scholar 

  49. Gao Z, Zhang J, Li L, Shen L, Li Q, Zou Y, Du X, Zhao Z. Heat shock proteins 27 and 70 contribute to the protection of Schisandrin B against d-galactosamine-induced liver injury in mice. Can J Physiol Pharmacol. 2016;94(4):373–8.

    Article  CAS  PubMed  Google Scholar 

  50. Ishaq M, Ojha R, Sharma K, Sharma G, Singh SK, Majumdar S. Functional inhibition of Hsp70 by Pifithrin-mu switches Gambogic acid induced caspase dependent cell death to caspase independent cell death in human bladder cancer cells. Biochim Biophys Acta. 2016;1863(11):2560–73.

    Article  CAS  PubMed  Google Scholar 

  51. Nosareva OL, Ryazantseva NV, Stepovaya EA, Shakhristova EV, Stepanova EA, Gulaya VS. The role of heat shock proteins 27 and 70 in redox-dependent regulation of apoptosis in Jurkat tumor cells. Biomed Khim. 2016;62(6):670–3.

    Article  CAS  PubMed  Google Scholar 

  52. Leiros GJ, Ceruti JM, Castellanos ML, Kusinsky AG, Balana ME. Androgens modify Wnt agonists/antagonists expression balance in dermal papilla cells preventing hair follicle stem cell differentiation in androgenetic alopecia. Mol Cell Endocrinol. 2017;439:26–34.

    Article  CAS  PubMed  Google Scholar 

  53. Magtanong L, Ko PJ, Dixon SJ. Emerging roles for lipids in non-apoptotic cell death. Cell Death Differ. 2016;23(7):1099–109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Rendl M, Lewis L, Fuchs E. Molecular dissection of mesenchymal-epithelial interactions in the hair follicle. PLoS Biol. 2005;3(11):1910–24.

    Article  CAS  Google Scholar 

  55. Schneider MR, Schmidt-Ullrich R, Paus R. The hair follicle as a dynamic Miniorgan. Curr Biol. 2009;19(3):R132–42.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by grants from National Natural Science Foundation of China (31772571, 31402038), China Postdoc Science Foundation (2014 M560810, 2015 T81059), as well as Key Research Program of Shaanxi Province (2017NY-072). Xiaolong Wang was supported by the Innovative Talents Promotion Plan in Shaanxi Province, and is a Tang Scholar at Northwest A&F University.

Availability of data and materials

All raw data have been submitted on NCBI, the accession numbers following as bellow, SRR6076964, SRR6076965, SRR6076962, SRR6076963, SRR6076960 and SRR6076961 for miRNA libraries; SRR6075308, SRR6075307, SRR6075306, SRR6075305, SRR6075310 and SRR6075309 for lncRNA libraries.

Author information

Authors and Affiliations

  1. College of Animal Science and Technology, Northwest A&F University, Yangling, 712100, China

    Guangxian Zhou, Danju Kang, Sen Ma, Xingtao Wang, Ye Gao, Yuxin Yang, Xiaolong Wang & Yulin Chen

Authors
  1. Guangxian Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Danju Kang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Sen Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xingtao Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Ye Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yuxin Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Xiaolong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yulin Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceived and designed the experiments: GX Z, YL C, XL W. Performed the experiments: GX Z, DJ K. Analyzed the data: GX Z. Contributed reagents and materials: S M, XT W. Contributed comments: Y G, XL W, YX Y. Wrote the paper: GX Z, DJ K, XL W. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Xiaolong Wang or Yulin Chen.

Ethics declarations

Ethics approval and consent to participate

All sampling procedures in this experiment were in accordance with approved guidelines of the Animal Care and Use Committee of the Northwest A&F University (Approval ID: 2014ZX08008–002).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Primers used in this study. (XLSX 10 kb)

Additional file 2:

Summary data output for each sample (DOCX 17 kb)

Additional file 3:

Figure S1. Raw reads classification for each sample (lncRNAs); Figure S2. Mapped region for each sample (lncRNAs);Figure S3. Novel lncRNAs filter. (a) Merged all sample transcripts and remove the chain direction uncertain transcript. Step1, the number of transcripts possessing no less than 2 exons; step 2, the number of transcripts possessing more than 200 bp length; step 3, the left number of transcripts getting rid of known lncRNAs; step 4, the left transcripts which FPKM ≥0.5; step 5, three coding potential screening tools (Coding-Non-Coding-Index, CNCI; Coding Potential Calculator, CPC and Pfam Scan, PFAM) applied to predict novel lncRNAs. (b) A venn diagram for step 5.;Figure S4. LncRNAs targeted genes GO analysis. (a–c) The significant molecular function, biological process and cellular component were enriched by lncRNAs targeted mRNAs in goat skins, and the DAG is the graphical display of GO enrichment results with candidate genes. (d) The number of genes in GO term were showed in histograph. Figure S5. GO analysis of lncRNA-miRNA target mRNAs. (a–c) The significant molecular function, biological process and cellular component were enriched by lncRNA-miRNA-mRNAs in goat skins, and the DAG is the graphical display of GO enrichment results with candidate targeted genes. (d) The number of genes in GO terms were showed in histograph. Figure S6. GO analysis of DE genes. (a–c) The significant molecular function, biological process and cellular component were enriched by DE mRNAs in goat skins, and the DAG is the graphical display of GO enrichment results with candidate targeted genes. (d) The number of genes in GO terms were showed in histograph. Figure S7. The top 20 KEGG pathways of hair cycle in skin. When the rich factor is greater, the Q-value is closer to zero, and the number of genes is greater, then the enrichment is more significant. (PDF 2925 kb)

Additional file 4:

FPKM for total transcripts and novel transcripts information. (XLSX 5602 kb)

Additional file 5:

TMP for total miRNAs and novel miRNAs information. (XLSX 108 kb)

Additional file 6:

The information for DGE lncRNAs. (XLSX 25 kb)

Additional file 7:

The information for DGE coding mRNAs. (XLSX 353 kb)

Additional file 8:

The information for DGE miRNAs. (XLSX 14 kb)

Additional file 9:

The enriched pathway terms. (XLSX 100 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, G., Kang, D., Ma, S. et al. Integrative analysis reveals ncRNA-mediated molecular regulatory network driving secondary hair follicle regression in cashmere goats. BMC Genomics 19, 222 (2018). https://doi.org/10.1186/s12864-018-4603-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4603-3

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us