Skip to main content

Integrated analysis of lncRNAs and mRNAs by RNA-Seq in secondary hair follicle development and cycling (anagen, catagen and telogen) of Jiangnan cashmere goat (Capra hircus)

Abstract

Background

Among the world’s finest natural fiber composites is derived from the secondary hair follicles (SHFs) of cashmere goats yield one of the world's best natural fibres. Their development and cycling are characterized by photoperiodism with diverse, well-orchestrated stimulatory and inhibitory signals. Long non-coding RNA (lncRNAs) and mRNAs play important roles in hair follicle (HF) development. However, not many studies have explored their specific functions in cashmere development and cycling. This study detected mRNAs and lncRNAs with their candidate genes and related pathways in SHF development and cycling of cashmere goat. We utilized RNA sequencing (RNA-Seq) and bioinformatics analysis on lncRNA and mRNA expressions in goat hair follicles to discover candidate genes and metabolic pathways that could affect development and cycling (anagen, catagen, and telogen).

Results

We identified 228 differentially expressed (DE) mRNAs and 256 DE lncRNA. For mRNAs, catagen and anagen had 16 upregulated and 35 downregulated DEGs, catagen and telogen had 18 upregulated and 9 downregulated DEGs and telogen and anagen had 52 upregulated and 98 downregulated DEGs. LncRNA witnessed 22 upregulated and 39 downregulated DEGs for catagen and anagen, 36 upregulated and 29 downregulated DEGs for catagen and telogen as well as 66 upregulated and 97 downregulated DEGs for telogen and anagen. Several key genes, including MSTRG.5451.2, MSTRG.45465.3, MSTRG.11609.2, CHST1, SH3BP4, CDKN1A, GAREM1, GSK-3β, DEFB103A KRTAP9–2, YAP1, S100A7A, FA2H, LOC102190037, LOC102179090, LOC102173866, KRT2, KRT39, FAM167A, FAT4 and EGFL6 were shown to be potentially important in hair follicle development and cycling. They were related to, WNT/β-catenin, mTORC1, ERK/MAPK, Hedgehog, TGFβ, NFkB/p38MAPK, caspase-1, and interleukin (IL)-1a signaling pathways.

Conclusion

This work adds to existing understanding of the regulation of HF development and cycling in cashmere goats via lncRNAs and mRNAs. It also serves as theoretical foundation for future SHF research in cashmere goats.

Peer Review reports

Background

In the 1970s, Xinjiang became China’s first international cashmere market, and it has since become a cashmere powerhouse. China provides a distinct husbandry resource called cashmere. It produces 12,000 tons of cashmere per year, accounting for more than 70% of global output [12]. Cashmere goat fleece have two coats; wool is made by primary hair follicles (PHF) and cashmere is produced by secondary hair follicles (SHF) [3,4,5]. Cashmere is a delicate (dehaired) downy undercoat shed by goats that lies beneath the outer thick hair shaft to provide an adaptability to cold [16]. Primary and secondary follicles differ particularly in terms of morphogenesis, cycle development, and positioning of hair follicles. The SHF does not have sebaceous gland like PHF [4] and follicle diameter and dermal papilla (DP) size are larger in PHF than the SHF [7]. Jiangnan cashmere goat is widely known for producing fine and silky hair fiber with great quality, economic value and high yield from SHF. Its extraordinary fiber features are related to its genetics and breeding qualities [8]. Cashmere is eight times warmer and just as soft as sheep wool, providing excellent insulation [9]. It remains one of the most valuable natural products in sales. While current output is dropping, global demand for cashmere in the textile and fabric business is increasing [1011]. Consequently, breeding high-quality and super-fineness cashmere has become a pressing issue in their breeding sector. SHF developmental characteristics influence cashmere yield and quality [12]. As a result, studies on the molecular mechanisms regulating cashmere goat HF development have become a research focus. As cashmere development is characterized by photoperiodic dependency, which occurs with cyclic regularity of SHF in quick growth (anagen), slow regression (catagen), and rest stage (telogen) [361314], exploring the main genes, signaling pathways and the complex regulatory mechanisms that enable their formation and development are important and must be investigated. Several key pathways involved in cashmere goat HF development and cycling include tumor necrosis factor (TNF), MAPK signaling pathway, WNT pathway, fibroblast growth factor (FGF) family, bone morphogenetic protein (BMP) family, transforming growth factor (TGF) family, Sonic hedgehog (SHH) conduction pathway and NOTCH signal transduction pathway [315,16,17,18]. However, not many signaling pathways that inhibit or enhance SHF development have been extensively investigated, and certain pathways are yet to be uncovered and comprehended. Current studies looking beyond protein-coding genes have concluded that non-coding RNA (ncRNA) such as microRNA (miRNA), natural antisense transcripts (NAT) and long non-coding RNA (lncRNA) may be more specific as biomarkers for certain applications than protein coding genes [13]. Among these ncRNAs, lncRNAs are abundant and also play crucial roles in cell development, cycle, proliferation and differentiation [18,19,20]. They perform functions, including transcriptional activation, protein coding gene silencing and mRNA or miRNA interaction to regulate their actions [2122]. While the understanding of the diverse regulatory mechanisms of ncRNAs is far from thorough, new studies have shown specific lncRNA pathways in HF morphogenesis in goats and sheeps [618]. A previous study in the expression profiles of lncRNAs and mRNAs throughout sheep fetal and postnatal hair follicle demonstrated that the integrative function of lncRNA and their target genes regulate HF development [23]. Another study used strand specific RNA sequencing (ssRNA-Seq) to ascertain the roles of lncRNAs and mRNAs in sheep skin at the on-set of SHF development and it was concluded that the key differentially expressed genes (DEGs) and lncRNAs may affect the molecular mechanisms involved in HF initiation [18]. Again, a study of the regulatory relationships of lncRNAs, miRNAs and mRNAs in the goat HF cycle was performed mainly in the catagen and anagen phases. It was revealed that lncRNAs and miRNAs directly correlate in HF cycling and catagen inducer factors transforming growth factor beta 1 (TGFβ1) and brain-derived neurotrophic factor (BDNF) were regulated by miR-873 and lnc108635596 [24]. A comprehensive analysis of DE mRNAs and DE lncRNAs between Inner Mongolia cashmere goats and Liaoning cashmere goats revealed that lncRNA XLOC_008679 and KRT35 target gene could be candidates for regulating cashmere fineness [5]. However, the roles of lncRNAs and mRNAs in cashmere development and cycling have not been thoroughly explored.

In this current study, the SHF of Jiangnan cashmere goat serve as an ideal model for investigating hair biology because of the circannual rhythmicity and synchronized growth of the fiber. We used RNA sequencing (RNA-seq) to investigate the profiles of lncRNAs and mRNAs in Jiangnan cashmere goat SHF in anagen, catagen and telogen stages as well as explore the potential interactions with regulatory networks of the related DEGs. This study will add more to our understanding of lncRNA and mRNAs in HF morphogenesis and contribute to the improvement of cashmere goat breeding. It will also serve as beneficial reference for Jiangnan cashmere goat genetic breeding as well as contribute to the annotation of the goat genome.

Materials and methods

Animal ethics

All animal experiments were carried out in strict accordance with the instructions provided by the Animal Care and Use Committee of Xinjiang Academy of Animal Science (Approval number 2020008). Effective procedures were implemented to reduce pain and distress; overall health, zoonotic infections, and pathogenic microbial infectious diseases were all thoroughly controlled and monitored [6, 25].

Experimental animal and sample collection

Experimental cashmere goats with high production features were collected from breeding farm in breeding center of Wenshu County, Aksu Prefecture, Xinjiang Province. The cashmere goats were fed in accordance with the standards of Xinjiang. Three female adults were chosen (two-year-old Jiangnan cashmere goat with a coefficient of relationship of < 0.125). Scapular skin tissues were collected from the goats in vivo. From each piece of skin tissue, approximately two cm2 and was immediately frozen by liquid nitrogen and stored at − 80 °C until further analysis [6].

Total RNA isolation, library preparation and sequencing

Total RNA was isolated and purified (Invitrogen, USA) and lncRNA and mRNA libraries were prepared as described before [6, 26]. The integrity of total RNA was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies Inc., USA) and sample with RNA Integrity number (RIN) values higher than 7.0 were used for sequencing. One microgram of RNA was used as input material for the RNA sample preparation. Subsequently, 9 cDNA libraries for the three developmental stages were constructed using rRNA-depleted RNA with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s instructions and each sample type contained three biological replicates. Libraries were sequenced using Illumina Hiseq 4000 platform (LC Bio, Hangzhou, Zhejiang, China) and 150-bp paired-end reads were generated. The raw data were assessed for quality using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Clean reads were assessed by removing adapters and poly-N > 10%, and low-quality reads in the raw data using the Fastx toolkit. Clean reads were mapped to goat reference sequences (Capra hircus ARS1.93) using Hisat [6]. Reference genome and annotation files were obtained from the Ensembl browser (https://www.ncbi.nlm.nih.gov/genome/?term=capra+hircus+ars1). Cufflinks software was used to assemble transcripts [27].

Coding potential and conserved analysis of lncRNAs

Based on the Cufflinks splicing results, lncRNA was chosen as final candidate for further analysis. The following parameters were used specifically; transcripts with transcript length ≥ 200 bp and exon number ≥ 2 were selected. The candidate transcripts were screened, using the class-code information from Cuffcompare. Cufflinks determined the transcript with coverage > 3 in at least one sample. Finally, we used coding potential analysis method, Coding-Non-Coding Index (CNCI), Coding Potential Calculator (CPC), PFAM and CPAT to choose the valid transcripts [6].

Analysis of differentially expressed (DE) lncRNAs and mRNAs

Gene abundance was expressed as fragment per kilobase of exons per million reads (FPKM) and lncRNAs and mRNAs in each sample was calculated using Stringtie software [28]. TMM algorithm was used for normalization. The R package edgeR [29] was used for differential expression analysis of mRNA and lncRNA. DE RNAs with p-values < 0.05 and |log2FoldChange| ≥ 1 served as thresholds to assess the statistical significance of mRNA and lncRNA expression differences.

Enrichment analysis and protein-protein interaction (PPI) of DE lncRNAs and mRNAs

GOseq R package was used for Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The Benjamini-Hochberg (BH) method [30] was used to adjust the significant p-values. GO terms and pathways with p-values < 0.05 were considered significantly enriched. The GO database was used to conduct functional annotations on DE transcripts, which were classified into three categories: biological processes (BPs), molecular function (MF), and cellular component (CC). Candidate genes encoded functional PPI network and were analyzed with STRING Genomics 11.0 database [31]. The network was visualized using the Cytoscape 3.6.0 program [32].

Correlation and co-expression analysis of lncRNAs and mRNAs

LncRNA target gene network was used to analyze the roles of lncRNAs. Cytoscape 3.6.0 software was used to create networks of lncRNA and target genes and they were predicted in cis and trans. The Pearson correlation coefficients (PCCs) between lncRNAs and their linked mRNAs were calculated to identify co-expression analysis. Certain indicators demonstrated strong correlation, and the important role of lncRNA or mRNA in the network is determined by the degree of correlation. The absolute value of Pearson’s correlation coefficient (r2) > 0.90 and p-value < 0.05 were set for subsequent networking of lncRNA and mRNA [6].

Quantitative reverse transcription PCR validation

Total RNAs were extracted from goat skin samples from the three groups and used for quantitative reverse transcription PCR (RT-qPCR). First strand cDNA was obtained using a One-step cDNA synthesis kit (Bio-Rad, USA) according to manufacturer’s instructions and was followed by subjection to quantification with GAPDH as endogenous control using standard SYBR Green PCR Kit (Bio-Rad) on the Bio-Rad CFX96 Touch™ Real Time PCR Detection System. The quantitative PCR was carried out following protocols: 95 °C for 35 s, followed by 42 cycles of 95 °C for 5 s and the optimized annealing temperature at 60 °C for 35 s and extension at 72 °C for 35 s. The primers used are listed in Table 12. 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 GAPDH. Student’s t-test were used for statistical analysis and p-value < 0.05 as significant. Values were expressed as means ±SD, * P < 0.05, ** P < 0.01.

Results

Overview of high-throughput sequencing

After extracting total RNA from three Jiangnan cashmere goats in the anagen (An), catagen (Cn), and telogen (Tn) phases, we constructed nine transcriptome libraries of cashmere goat skin samples and sequenced the RNA. Quality control findings revealed a total of 109.40Gb clean data and the percentage of Q20 of each sample was > 98.18% and the percentage of Q30 base of each sample was > 94.62%. The GC content ranged between 46.09 and 47.42% (Table 1). The numbers 1, 2 and 3 in the sample names signify the 1–3 cashmere goats. According to the quality control results, the sequencing results were reliable and adequate for further data processing. Based on the comparative results, differential splicing predictions assessment, gene structure optimization analysis, and novel gene extraction were conducted to uncover 10,054 new genes, of which 2587 were functionally annotated. 228 DEGs were discovered based on the level of gene expression in various samples and their functional annotation and enrichment analysis were performed. There were 27,740 lncRNAs discovered, with 256 DE lncRNAs.

Table 1 RNA-Seq quality control result

After mapping the clean reads to goat reference genome (Capri hircus ARS1) (Table 2), comparison of each growth cycles proceeded. The three comparing groups were: (1) anagen vs catagen, (2) catagen vs telogen and (3) telogen vs anagen. It was demonstrated by limiting the q-value to < 0.05. Most of the DEGs were downregulated in the first and third group (Table 3). Functional annotations of DEGs and statistical results of the genes annotated by the gene for each differential expression is shown in Table 4.

Table 2 Reference genome alignment read statistics
Table 3 Statistical table of the number of DEGs
Table 4 Functional annotation for DEGs among the different stages

The DEGs were depicted in a Venn diagram (Fig. 1). The number of DEGs was highest in telogen and anagen, and gene expression varied significantly in the telogen stage compared to the other phases.

Fig. 1
figure 1

Identification of mRNAs in skin tissues of Jiangnan cashmere goat. The number of mRNAs that are unique to each comparison group, as well as the number of mRNAs that are shared by all comparison groups. Each circle indicates a unique combination of differential expression

The mRNAs were grouped using hierarchical clustering based on their similarities in gene expression patterns across the three groups (Fig. 2). Gene expression pattern in the telogen phase appeared to be significantly different compared to the other phases.

Fig. 2
figure 2

Cluster heatmap of DE mRNA on the basis of their expression values. Orange and green indicate high and low expressions, respectively. The colors denote differential expression levels (log2 (fold change) ≥1 and p-value< 0.05)

Volcano maps were plotted to analyze the differences in gene expression levels between groups as well as statistical significance of the differences Fig. 3.

Fig. 3
figure 3

Volcano plots of DE mRNAs at different cycles of SHF development. a catagen vs anagen (b) catagen vs telogen (c) telogen vs anagen. The y-axis indicates the −log10(FPKM+ 1) values, the green points in the figure represent the down-regulated, the red dot represents the up-regulated, and the black represents the non-difference expression

Enrichment of DEGs

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) analyses were performed to better understand the biological processes and pathways involved in cashmere HF development. Following the comparison of DEGs across groups, the top enriched GO terms (p-value < 0.05) were analyzed. Among the groups the DEGs were related to biological processes like fatty acid elongation, polyunsaturated fatty acid (GO:0034626), long-chain fatty acid catabolic process (GO:0042758), fatty acid elongation, saturated fatty acid (GO:0019367), long-chain fatty-acyl-CoA biosynthetic process (GO:0035338) and unsaturated fatty acid biosynthetic process (GO:0006636), keratinocyte activation (GO:0032980), medium-chain fatty acid catabolic process (GO:0051793), skin development (GO:0043588), pointed-end actin filament capping (GO:0051694), collagen catabolic process (GO:0030574) and molecular function terms like calcium ion binding (GO:0005509), transition metal ion binding (GO:0046914), fatty acid elongase activity (GO:0009922), high affinity glutamate transmembrane transport activity (GO:0005314) as well as the cellular function terms like lipid droplet (GO:0005811), keratin filament (GO:0045095), suggesting their developmental and maintenance functions in cashmere HF cycle. The top 20 GO terms among the various cycles are shown in (Fig. 4).

Fig. 4
figure 4

Enriched GO terms for DE mRNAs (p-value < 0.01). a catagen vs anagen (b) catagen vs telogen (c) telogen vs anagen

According to the KEGG pathway, there were significant changes in the physiological processes in the various phases of the SHF development and cycling. Catagen vs Anagen showed enrichment of 35 KEGG pathways, including Arachidonic acid metabolism, glycolipid metabolism, fatty acid elongation, fatty acid metabolism, folate biosynthesis and Hippo signaling pathway. Catagen vs telogen showed enrichment of 15 KEGG pathways, including folate biosynthesis, mineral absorption, IL-17 signaling pathway, glutamatergic synapse and arachidonic acid metabolism. Telogen vs anagen showed enrichment of 143 KEGG pathways, including Hippo signaling pathway, WNT signaling pathway, JAK-STAT signaling pathway, NF-kappa B signaling pathway, Th17 cell differentiation, Th1 and Th2 cell differentiation, PI3K-Akt signaling pathway, VEGF signaling pathway, Arachidonic acid metabolism and IL-17 signaling pathway. The top 20 KEGG pathways [33,34,35] among the different cycles were shown in (Fig. 5).

Fig. 5
figure 5

Enriched KEGG pathway for DE mRNAs (p-value < 0.05). a catagen vs anagen (b) catagen vs telogen (c) telogen vs anagen. The ratio of mRNAs enriched in the pathway among those identified in the pathway is indicated by the rich factor

Identification of lncRNAs

37,386 lncRNAs were discovered and their genomic contexts were examined based on their similarities to local protein-coding genes. They were classified into (13749) long intergenic RNAs (lincRNAs), (2469) Antisense-lncRNA, (11103) Intronic-lncRNA, (419) sense lncRNA (Fig. 6a). After excluding candidate lncRNAs with coding potential using the software CNCI, CPC, CPAT and PFAM-scan protein domain analysis, 27,740 lncRNAs were identified (Fig. 6b).

Fig. 6
figure 6

Identification of long noncoding RNAs (lncRNAs) in Jiangnan cashmere goat skin tissues. a Bar chart showing the four categories of lncRNA. b Venn diagram showing the number of lncRNAs with coding potential analysis by CNCI, CPC, CPAT and PFAM

The number of lncRNA expressed by the differences in each group is demonstrated in Table 5. Tn-1,Tn-2,Tn-3 vs An-1,An-2,An-3 had the highest DEGs.

Table 5 Differentially expressed lncRNA

The number of lncRNAs with unique and common differences between the comparison groups were plotted. Each set of DE lncRNA was demonstrated in Venn diagram (Fig. 7).

Fig. 7
figure 7

Venn diagram showing the number of DE lncRNA with differences and similarities among the comparison groups. A catagen vs anagen B catagen vs telogen C telogen vs anagen

The DE lncRNAs were grouped using hierarchical clustering based on their similarities in gene expression patterns across the three groups (Fig. 8).

Fig. 8
figure 8

Clustered heat map of DE lncRNAs on the basis of their expression values. Orange and green indicate high and low expressions, respectively. The colors denote differential expression levels (log2 (fold change) ≥1 and p-value< 0.05)

DE lncRNA functional classification analysis in Jiangnan cashmere goats

EdgeR, GO and KEGG analyses were used to investigate the functional clustering characteristics of DE lncRNAs between cashmere goat seasons. The lncRNAs cis target genes were enriched in important GO terms related to HF development. Among the groups the lncRNAs cis target genes were related to biological processes such as cornification (GO:0070268), epidermal cell differentiation (GO:0009913), keratinization (GO:0031424), lipid biosynthesis process (GO:0008610), keratinocyte differentiation (GO:0030216), regulation of BMP signaling pathway (GO:0030510), lipid modification (GO:0032102), epidermal growth factor receptor signaling pathway (GO:0007173), positive regulation of ERK1/2 cascade, Wnt signaling pathway (GO:0016055), integrin-mediated signaling pathway, Notch signaling pathway (GO:0007219), positive regulation of vascular endothelial growth factor production (GO:0010575) and hippo signaling (GO:0035329). The cellular component terms showed enrichment in keratin filament (GO:0045095), endoplasmic reticulum membrane (GO:0005789), integrin complex (GO:0008305), growth cone (GO:0030426) and actin filament (GO:0005884). The molecular functions among the various cycles showed enrichment in transforming growth factor beta binding (GO:0050431), growth factor receptor binding (GO:0070851), NF-kappaB binding (GO:0051059), collagen binding (GO:0005518), beta-catenin binding (GO:0008013), actin filament binding (GO:0051015) and microtubule binding; indicating their potential developmental functions and maintenance activities in cashmere cycling. The top 20 GO terms among the various cycles are shown in (Fig. 9).

Fig. 9
figure 9

Enriched GO terms for lncRNAs cis-target genes (p-value < 0.01) (a) catagen vs anagen (b) catagen vs telogen (c) telogen vs anagen

According to the KEGG pathway, there were significant changes in the physiological processes occurring at the various phases of cashmere goat hair follicle development and cycling. Some of the observed KEGG pathways belonged to conventional techniques linked with the cashmere HF cycle. Catagen vs Anagen showed enrichment of 108 KEGG pathways including pl3k-Akt signaling pathway, Hippo signaling pathway, IL-17 signaling pathway, fatty acid metabolism, Notch signaling pathway, TGF-beta signaling pathway, glycerolipid metabolism, fat digestion and absorption, MAPK signaling pathway, glutamatergic synapse and cell cycle. Catagen vs telogen showed enrichment of 171 KEGG pathways including MAPK signaling pathway, Th1 and Th2 cell differentiation, PPAR signaling pathway, Notch signaling pathway, NF-Kappa B signaling pathway, PI3K-Akt signaling pathway, Wnt signaling pathway, IL-17 signaling pathway and Hedgehog signaling pathway. Telogen vs anagen showed enrichment of 222 KEGG pathways including JAK-STAT signaling pathway, Th17 cell differentiation, Hedgehog signaling pathway, Arachidonic acid metabolism, NF-Kappa B signaling pathway, Notch signalling pathway, pl3k-Akt signaling pathway, IL-17 signaling pathway, Hippo signaling pathway and Wnt signaling pathway. The top 20 KEGG pathways [33,34,35] with the highest representation of lncRNAs cis-target genes are shown in (Fig. 10).

Fig. 10
figure 10

Top 20 significantly (p-value < 0.05) enriched KEGG pathway for lncRNAs cis-target genes (a) between catagen and anagen (b) between catagen and telogen (c) between telogen and anagen. Rich factor represents the ratio of DE lncRNA cis-target genes enriched in the pathway among genes annotated in the pathway

The lncRNAs trans target genes were enriched in important GO terms relating to HF development. Among the groups the lncRNAs trans target genes are related to biological processes such as positive regulation of NF-kappaB transcription factor activity, cell cycle arrest (GO:0007050), positive regulation of B cell differentiation (GO:0045579), fatty acid beta-oxidation (GO:0006635), regulation of ERK1/2 cascade (GO:0070372), lipid transport (GO:0006869), positive regulation of NF-KappaB transcription factor activity (GO:0051092), lipid transport (GO:0006869), fatty acid transport (GO:0015908), positive regulation of MAPK cascade (GO:0043410), positive regulation of canonical Wnt signaling pathway, germ cell development (GO:0007281) and integrin mediated signaling pathway (GO:0007229). The cellular component terms showed enrichment in autophagosome (GO:0005776), keratin filament (GO:0045095), integrin alphav-beta8 complex (GO:0034686), actin cytoskeleton (GO:0015629), cilium (GO:0005929), microtubule cytoskeleton (GO:0015630) and sarcolemma (GO:0042383). The molecular functions among the various cycles showed enrichment in magnesium ion binding (GO:0000287), protein serine/threonine kinase activity (GO:0004674), transcription factor activity, sequence-specific DNA binding (GO:0003700), transferase activity, transferring acyl groups other than amino-acyl groups (GO:0016747), ATP binding (GO:0005524), calcium ion binding (GO:0005509), magnesium ion transmembrane transporter activity (GO:0015095), zinc ion binding (GO:0008270), metal ion binding (GO:0046872) and protein binding (GO:00055115). The top 20 GO terms among the various cycles are shown in (Fig. 11).

Fig. 11
figure 11

Enriched GO terms for lncRNAs trans-target genes (p-value < 0.01) (a) catagen and anagen (b) catagen and telogen (c) telogen and anagen

According to the KEGG pathways, there were significant changes in the physiological processes occurring at various phases of cashmere goat hair follicle development and cycling. Some of the observed KEGG pathways belonged to conventional pathways associated with cashmere hair follicle cycle. Catagen vs Anagen showed enrichment of 118 KEGG pathways including Pl3K-Akt signaling pathway, MAPK signaling pathway, WNT signaling pathway, Th17 cell differentiation, Th1 and Th2 cell differentiation and glycerolipid metabolism. Catagen vs telogen showed enrichment of 93 KEGG pathways including Pl3K-Akt signaling pathway, Th17 cell differentiation, Th1 and Th2 cell differentiation, WNT signaling pathway, MAPK signaling pathway and Glutathione metabolism. Telogen vs anagen showed enrichment of 139 KEGG pathways including WNT signaling pathway, MAPK signaling pathway, NF-kappaB signaling pathway, Th17 cell differentiation, Th1 and Th2 cell differentiation, cell cycle, Arachidonic acid metabolism, Glutamatergic synapse, Pl3K-Akt signaling pathway and calcium signaling pathway. The top 20 KEGG pathways [33,34,35] with the highest representation of lncRNAs trans-target genes are shown in (Fig. 12).

Fig. 12
figure 12

Top 20 significantly (p-value < 0.05) enriched KEGG pathway for lncRNAs trans-target genes (a) between catagen and anagen (b) between catagen and telogen (c) between telogen and anagen. Rich factor represents the ratio of DE lncRNA trans-target genes enriched in the pathway among genes annotated in the pathway

DE lncRNA screening was carried out and the top ten upregulated and downregulated genes, as well as their potential target genes were chosen in Cn1, Cn2, Cn3 vs An3, An2, An1 (Table 6), Cn1, Cn2, Cn3 vs Tn1, Tn2, Tn3 (Table 7) and Tn1,Tn2,Tn3 vs An1,An2,An3 (Table 8). The DE lncRNA’s targets were predicted using both cis and trans principles. The Fold Change ≥2 and FDR < 0.05 served as screening criteria. The fold change multiple represents the ratio of expression between two samples (groups). We found that some lncRNAs, MSTRG.46401.1, MSTRG.4712.2, MSTRG.48192.1 predominantly have trans target genes in the KRT family and is directly associated to the growth of goat hair follicles, and KRT dysregulation could result in hair diseases Table 6. MSTRG.48150.1, MSTRG.48153.1, MSTRG.48154, MSTRG.48155 have cis target genes mainly WNT8B, which is a member of WNT gene family, involved in regulating hair follicles morphogenesis.

Table 6 DE lncRNA of Cn1, Cn2, Cn3 vs An3, An2, An1
Table 7 DE lncRNA of Cn1, Cn2, Cn3 vs Tn1, Tn2, Tn3
Table 8 DE lncRNA of Tn1,Tn2,Tn3 vs An1,An2,An3

Characteristic analysis of lncRNA and mRNA

Transcription length, exon number, ORF length, expression level, variable shear isomer, and interactive analysis were used to analyze the differences and characteristics of lncRNA and mRNA in terms of understanding the structural features of genes in cashmere goats (Fig. 13). The length distribution study revealed that lncRNAs were longer than mRNAs. The prediction of cis and trans target genes for DE lncRNAs were used to provide better understanding to the functions of the lncRNAs at various stages of development.

Fig. 13
figure 13

Comparison of genomic primary differences and expression levels between lncRNAs and mRNAs in cashmere goat hair follicle. a The expression levels were indicated by log10(FPKM+ 1) in the lncRNAs and mRNAs. b Comparison of variable shear isomers for lncRNAs and mRNAs

The mRNAs compared with lncRNAs expression levels and mRNAs and lncRNAs variable shear isomer comparative analysis are shown in Fig. 13.

To understand the differences in architecture and expression levelsbetween the mRNAs and lncRNAs, different analyses were done. The mRNAs compared with lncRNAs length, lncRNAs exons and ORF length of lncRNAs Fig. 14.

Fig. 14
figure 14

Comparison of genomic architecture and expression level between lncRNAs and mRNAs in cashmere goat hair follicle. a Distribution of transcript length in lncRNAs and mRNAs in goat skin. The horizontal axis indicates the length of transcripts and the vertical axis represents density (b) Distribution of the number of exons in the mRNAs and lncRNAs. c Distribution of the number of open reading frames (ORFs) in the mRNAs and lncRNAs. The ORFs were identified using Estcan

DE lncRNA and DE mRNA interactive analysis

The expression of lncRNAs and mRNAs in different interactive analysis, in relation to the similarities and differences between the RNAs in different cycles.

Volcano maps were plotted to analyze the relationships in gene expression levels between groups as well as statistical significance of those relationships. The relationships between the groups were expressed in the volcano maps in Fig. 15.

Fig. 15
figure 15

Volcano plot of differentially expressed (DE) lncRNA and mRNA (log2(fold change) ≥1 and p-value < 0.05). A Catagen vs Anagen (B) Catagen vs Telogen (C) Telogen vs Anagen. The y-axis indicates the −log10(FPKM+ 1) values. Red represent upregulated lncRNA, green represents downregulated lncRNA, orange represents upregulated gene, blue represents downregulated gene, and black represents non-differential expression gene

MA interactive diagrams were plotted to analyze the relationships in gene expression levels between groups. The relationships between the groups are expressed in the MA diagram in Fig. 16.

Fig. 16
figure 16

MA interactive diagram of differences in the expression of lncRNA and mRNA (A) Catagen vs Anagen (B) Catagen vs Telogen (C) Telogen vs Anagen. The x-axis indicates log2 (FPKM). Red represents upregulated lncRNA, green represents downregulated lncRNA, orange represents upregulated genes, blue represents downregulated genes, and black points represent genes that express no significant difference

The sequencing results were validated using RT-qPCR. Figure 17 depicts the end results. It was normal that the expression of LOC1021179881 and LOC102181431 between RNA-Seq and qPCR were inconsistent. Some genes/ lncRNAs were low-abundant while sensitivity of RNA-Seq and qPCR were different which caused the difference in expression.

Fig. 17
figure 17

QPCR validation of RNA-Seq. Eight mRNAs and lncRNAs were selected for validation including LOC1021179881, LOC102181431, LOC100861181, MSTRG.55192.23, MSTRG.31147.1, KRTAP11–1, LUM and GLYATL2. The x-axis indicates the three main cycles of SHF and the y-axis indicates lncRNA and mRNA relative expression level (mean ± SEM)

Discussion

Cashmere goats exist as a type of cyclical molting animals used to study photoperiod-dependent hair follicle morphological changes [9]. Their SHFs account for the important part of the skin that produce fine variable cashmere fleece in repetitive cycles, traversing the phases of growth (anagen), regression induced by apoptosis (catagen) and relative quiescence (telogen) [3924]. The recurrent cycles are regulated by extensive and well-orchestrated interactions between follicular epithelial cells and mesenchymal cells of dermal papilla (DP) to necessitate the dynamic synthesis of various stimulatory and inhibitory signals [3637]. WNT, Notch, BMP/TGF, NF-kappa B, HGF, Shh, and EDA/EDAR signaling pathways are among those implicated in critical signaling pathways which intricate interplay and interdependence of HF development and cycling [3238,39,40,41]. At least, three key developmental signaling pathways (WNT, SHH, and NF-kB/EDAR) are essential for hair follicle growth and maintenance [317]. LncRNAs do not encode protein [1920]. They are at least 200 bp long and include conserved secondary structures [38]. Numerous lncRNAs have been shown to have critical roles in biological processes including apoptosis, proliferation and differentiation [3842]. They interact with DNAs, RNAs, and proteins to regulate gene expression through different interconnected mechanisms including cell cycle regulation, splicing regulation, gene imprinting, mRNA degradation and stabilization, and translational regulation [22384344]. The general characteristics of lncRNAs and mRNAs reveal that, they are relatively comparable in terms of being 5 capped, spliced and polyadenylated and lncRNAs are much less evolutionarily conserved, contain fewer exons and are less abundantly expressed [38]. The specific function of lncRNAs and mRNAs in HF development and cycling have not been extensively analyzed and as a result remain poorly understood. Since RNA expression levels change overtime and across tissues, transcriptome sequencing is regarded ideal for exploring gene expression variations [45]. In this study, we used RNA-Seq technology to investigate the SHF development cycle in cashmere goats. We analyzed lncRNAs and mRNAs in SHF development among the three main cycles (anagen, catagen and telogen). We identified 256 DE lncRNAs, some of which were trans and cis target genes and 228 DE mRNA involved in SHF development and cycling. The main limitation of the study is that, goat genomic annotation is incomplete.

The randomly selected lncRNAs and mRNAs used for validating the RNA-Seq results revealed the following findings (Table 9). GLYATL2 is involved in fat synthesis [46] and linked to hair follicle morphogenesis and cycling as it is involved in lipid modification of signaling proteins such as Hedgehogs (HHs) and WNTs [47]. Lum regulate development, growth and healing through the functional activities of cytokines. Lum together with Lgals1 and Apoa1 collectively enhance the adult fibroblast competence for HF neogenesis [48]. Lum/ Lgals1/ Apoa1 may help promote Wnt signaling [48]. Lum serves as a key node in the regulatory network and is required for collagen fibrillogenesis in skin homeostasis [10]. LOC102181431 (KRTAP13–1) is part of PMG protein family which is characterized by anagen specific KAPs and play important role in phosphorylation which regulates the formation of keratin filament cell cycle [49]. KRT and KRTAP constituent the major structural proteins of the hair fiber and sheath and their products serve as essential commodity for fleece quality and cashmere goat hair morphology [5]. KRTAP gene expression modulate the fineness and density of cashmere as well as other certain features [50]. LOC102179881 (KRT33A) is primarily expressed in hair follicle and highly expressed in its cortex, where it interacts with keratinization and developmental biology pathways [51]. KRTAP11–1 genetic variation may regulate its expression, protein structure and posttranslational modifications and subsequently affect wool fiber structure and wool traits [52]. MSTRG.31147.1 has target gene PADI1 is involved in hair follicle differentiation and may take part in transcriptional regulation for hair follicle keratin synthesis [5354]. MSTRG.55192.23 targets COL1A1 is involved in processes and pathways that promote cyclic growth of HF [32].

Table 9 Primers used in RT-qPCR analysis

Functional annotation of DEGs were analyzed using KEGG and GO database. DEGs significantly enriched in GO analysis included HSD3B5, GTF2IRD1, GDAP10, GYS3, GSN, CDK2AP1, CD151, ADGRE5 etc. (Fig. 4). KEGG pathway analysis (Fig. 5) showed that HOXA4, LTα, COQ7, ASSP12, KCND2, FSHR, SF3B4, GTL11 were enriched in Hippo (ko04391), VEGF (ko04370), Pl3k-Akt (ko04151), NF-kappa B (ko04064), Wnt (ko04310) and Hedgehog (ko04340) signaling pathway [33,34,35]. These findings imply that the mRNAs may play key roles in the signaling networks across distinct hair follicle developmental stages of cashmere goats. COQ7 has been shown to regulate a number of developmental and physiological mechanisms such as apoptosis, cell cycle, gene expression, rhythmic behavior and aging [55]. It was found in NF-kappa B (ko04064) signaling pathway, exclusively in telogen and anagen phase and may function as a key developmental signaling pathways. HOXA4 and LTα were related to Hippo (ko04391)/ (ko04392) signaling pathway, respectively. LTα may activate NF-kappaB-signaling pathway and implicate the control of normal hair follicle development and the pathophysiology of certain hair disorders [5657].

In the Cn1, Cn2, Cn3 vs An3, An2, An1 phase, a variety of DEGs in hair follicle development were identified (Table 10). The results showed that epidermal differentiation complex (EDC) genes were upregulated. EDCs promote keratinocyte differentiation singly and/or combination, leading to epidermal and hair follicle differentiation. They include S100 family; LOC102191570, LOC108633164, S100A7A and function as innate immune regulatory elements [58]. S100A7A is mainly found in cytoplasm and nucleus of different types of cells involved in continuous regulation of cellular migration, invasion, differentiation and cell cycle progression [59]. The expression levels of S100A7A in normal skin of humans are usually low and in altered conditions, their expression increase [58]. The strong upregulation of S100A7A could be characterized by increased immune regulatory functions as a result of psoriasis and other hyperproliferative skin conditions, featuring impaired epidermal differentiation. The main regulatory pathways of S100A7A expression are NFkB/p38MAPK, caspase-1, and interleukin (IL)-1a [59]. GLYATL2 is related to fat synthesis and it was significantly upregulated. KRT4 upregulation suggests significant increase in cell proliferation to promote SHF development and also increase cashmere yield [60]. Elevated expression of SLC1A2 suggests that it may promote hair growth and keratinocyte proliferation [61]. Dysregulation of LOC102190037/LOC102179090/LOC102173866 may impair the essential role of lipid metabolism in goat skin functions, probably because lack of essential fatty acids could reduce water repulsion and increase trans-epidermal water loss [62]. They may exert their functions through signaling proteins such a Hedgehogs (HHS) and WNTS, which are essential for hair follicle development and cycling [47]. Downregulation of FA2H may disrupt synthesis of 2-hydroxylated glucosylceramide and wax diesters in sebaceous gland and cause delayed hair fiber from follicles, hyperproliferation and cyclic alopecia [63]. In Cn1, Cn2, Cn3 vs Tn1, Tn2, Tn3, AQP8 and KRT2 were highly functional and could play key roles in hair follicle development (Table 11). AQP8 play significant roles in transportation of water, urea, ammonia, hydrogen peroxide and gland secretion in hair follicle development [64]. Keratins are well-known for their important involvement in hair follicle cell proliferation and differentiation [4965] and are predominantly found in hair cortex [66]. KRT2 mainly participates in epidermis keratinization and its expression was lower during telogen and this causes inactive hair growth [24]. In Tn1,Tn2,Tn3 vs An1,An2,An3 a number of mRNAs played key roles in hair follicle development (Table 12). KRT39 was upregulated and have shown to be potentially significant for regulating hair follicle fiber growth and fineness [6]. Decreased expression of β-defensins have been associated with human atopic dermatitis. LOC102188015 (DEFB103A) was upregulated and may play significant roles in innate immunity against bacterial skin infections [67]. FAM167A is transcribed in many cell lines and tissues and is suggested to exert functions in immune pathways and cell motility. It may be involved in autoimmune pathogenesis of lupus which scar hair follicles, causing hair loss [68]. EGFL6 is mostly expressed from the epidermis and it has shown to be essential for proper alignment, touch response and av. integrin-enrichment of lanceolate complexes [69]. The downregulation of EGFL6 may affect stable hair follicle-lanceolate complex interface responsible for tactile sensation [70]. FAT4 has been demonstrated in human cancer, including melanoma. It performs significant function in suppressing tumor growth through activation of Hippo signaling and its genetic perturbations in hair follicle has been implicated in their development and cycling [71]. In this study, downregulation of FAT4 may induce cell cycle arrest [72].

Table 10 The top 10 upregulated and downregulated DEGs of Cn1, Cn2, Cn3 vs An3, An2, An1
Table 11 The top 10 upregulated and downregulated DEGs of Cn1, Cn2, Cn3 vs Tn1, Tn2, Tn3
Table 12 The top 10 upregulated and downregulated DEGs of Tn1,Tn2,Tn3 vs An1,An2,An3

LncRNA is known to play important roles in regulating development of cashmere goats and therefore serve as potential candidate genes for studying the molecular mechanisms of hair follicle morphogenesis. Despite the fact that they do not encode proteins, LncRNAs have the ability to regulate gene expression and stability of target genes [42]. The regulatory pattern characteristics of lncRNAs indicate that they may be classified as cis and trans. The cis regulate adjacent genes whereas trans regulation occur throughout the cell [20]. MSTRG.11609.2 target gene CHST11 has demonstrated to be positively regulated by transforming growth factor-beta (TGFβ) signaling pathway [73] and plays key roles in HF development. This study suggests that CHST11 promotes proliferation of dermal papillae for hair growth partly through activation of WNT/β-catenin signaling of HF cells. MSTRG.5451.2 target gene SH3BP4 is involved in indirect regulation of cell growth, proliferation and autophagy. It negatively regulates GTPase complex and amino acid dependent mTORC1 signaling [74]. The upregulation of mTORC1 signaling by the target gene also affect WNT/ β-catenin signaling and may cause loss of HF maintenance and impaired hair growth. Although, overactive WNT signaling leads to HF degeneration in anagen initiation [75]. MSTRG.41648.1 target gene FBLN5 binds integrins to modulate differentiation in the hair shaft and HF development [76]. MSTRG.52443.14 target gene CDKN1A and PNPLA4. CDKIs are potent inhibitors of cell cycle progression that are transcriptionally up-regulated in the bulge at all phases or at the start of HF quiescence [77]. CDKN1A may play important role in antagonizing the progression of HF development to cause cell cycle arrest at the catagen phase. PNPLA4 may influence epidermal homeostasis via transcription factor ligands such as retinoid and triglyceride [78]. MSTRG.45465.3 target gene GAREM1 plays important role in the regulation of cell proliferation. The GAREM1 is functionally important in ERK/MAPK signaling pathway and it is implicated in Y-27632 to induce proliferation of different stages ligament stem cells [7980]. It also regulates keratinocyte proliferation and differentiation [79]. MSTRG.840.3 target gene GSK3β and LOC108638286. The downregulation of GSK3β in this study, may be as a result of activation of WNT signaling. GSK-3β downregulation may inhibit the formation of the axin complex and indirectly elevates the level of β -catenin. Its inhibitors increase β-catenin levels in HF morphogenesis and SC differentiation in the hair cycle, hence initiating the transcription of a number of downstream genes [8182]. It has been shown that plucked human hair contained cells with high GSK-3β levels [81]. The downregulation of LOC108638286 (KRTAP9–2) may affect keratinization, as well as HF development and cycling [83]. Keratin-associated proteins (KRTAPs) aid in the development of rigid hair shaft by forming large disulfide bonds with keratin intermediate filament proteins [50]. A recent study revealed (KRTAP9–2) expression in all phases of the HF cycle. MSTRG.23343.1 target gene MCC is involved in the negative regulation of cell cycle progression. It serves as a candidate colorectal tumor suppressor gene and suppresses cell proliferation and WNT/β-catenin pathway in colorectal cancer cells [84]. Its role in HF development and growth needs further studies. MSTRG.32111.1 target BIRC3 and YAP1. BIRC3 regulates apoptosis, cell proliferation, cell invasion and metastasis. YAP1 is a transcriptional coactivator that regulates Hippo signaling pathway to control homeostasis, proliferation and differentiation [85]. It has shown to regulate SC activation and development cycling as well as regeneration of HFs in adult mice [85]. A previous study in mouse showed that K14 regulates early differentiation and proliferation through interactions with 14–3-3 adaptor proteins and YAP1. More studies should be done in the future to identify the specific interactions among certain genes identified in this study. Protein interaction network for DE lncRNA target genes demonstrated that ENSBTAG00000025441 (HSPA1A) and ENSBTAG00000004912 (ANKRD17) were key regulators of the cell cycle and innate immune defense.

Conclusion

In this study, we used RNA-Seq technique to screen and analyze the roles of DE mRNA and lncRNA among different cycles (anagen, catagen and telogen) in Jiangnan cashmere goat SHF. We identified DEGs in SHF potentially related to hair development and cycling. The DEGs included CHST11, SH3BP4, FBLN5, CDKN1A, PNPLA4, GAREM1, GSK3β, KRTAP9–2, BIRC3, GLYATL2, KRT4, SLC1A2, FA2H, KRT2, KRT39, FAM167A and EGFL6 and they were related to cell cycle, cell differentiation, proliferation, apoptosis, immune response, keratinization, proliferation, development, transcription regulation and cell growth. The results of this study could be used to facilitate further studies of the detailed molecular mechanisms in hair development and also serve as a foundation for future studies to potentially manipulate cashmere traits.

Availability of data and materials

The datasets presented in this study can be found in online repositories. The name of the repository/repositories and accession numbers can be found at Sequence Read Archive (SRA) repository, PRJNA778726.

References

  1. Jiang H-Z, et al. Industry status of Chinese cashmere goat and analysis of its prospects. Anim Husb Feed Sci. 2009;10:100–3.

    Google Scholar 

  2. Jin M, et al. Genetic signatures of selection for cashmere traits in Chinese goats. Animals (Basel). 2020;10(10):1905.

    Article  Google Scholar 

  3. Geng R, Yuan C, Chen Y. Exploring differentially expressed genes by RNA-Seq in cashmere goat (Capra hircus) skin during hair follicle development and cycling. PLoS One. 2013;8(4):e62704.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  4. Qiao X, et al. Discovery of differentially expressed genes in cashmere goat (Capra hircus) hair follicles by RNA sequencing. Genet Mol Res. 2016;15(3). https://doi.org/10.4238/gmr.15038589.

  5. Zheng YY, et al. An integrated analysis of cashmere fineness lncRNAs in cashmere goats. Genes (Basel). 2019;10(4):266.

    CAS  Article  Google Scholar 

  6. Fu X, et al. Integrated analysis of lncRNA and mRNA reveals novel insights into cashmere fineness in Tibetan cashmere goats. PeerJ. 2020;8:e10217–7.

    PubMed  PubMed Central  Article  Google Scholar 

  7. Zhu B, et al. Transcriptome sequencing reveals differences between primary and secondary hair follicle-derived dermal papilla cells of the cashmere goat (Capra hircus). PLoS One. 2013;8(9):e76282–2.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. Liu S, et al. Estimates of linkage disequilibrium and effective population sizes in Chinese merino (Xinjiang type) sheep by genome-wide SNPs. Genes Genomics. 2017;39(7):733–45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. Xing W, et al. Automatic identification of cashmere and wool fibers based on the morphological features analysis. Micron. 2020;128:102768.

    PubMed  Article  Google Scholar 

  10. Ge W, et al. Melatonin promotes cashmere goat (Capra hircus) secondary hair follicle growth: a view from integrated analysis of long non-coding and coding RNAs. Cell Cycle. 2018;17(10):1255–67.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. Nocelli C, et al. Shedding light on cashmere goat hair follicle biology: from morphology analyses to transcriptomic landascape. BMC Genomics. 2020;21(1):458.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. Nixon A, et al. Seasonal hair follicle activity and fibre growth in some New Zealand cashmere-bearing goats (Caprus hircus). J Zool. 1991;224(4):589–98.

    Article  Google Scholar 

  13. Wang S, et al. Integrated analysis of coding genes and non-coding RNAs during hair follicle cycle of cashmere goat (Capra hircus). BMC Genomics. 2017;18(1):767.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. Su R, et al. Screening the key genes of hair follicle growth cycle in inner Mongolian cashmere goat based on RNA sequencing. Arch Anim Breed. 2020;63(1):155–64.

    PubMed  PubMed Central  Article  Google Scholar 

  15. Zhang Y, et al. Comparative study on seasonal hair follicle cycling by analysis of the transcriptomes from cashmere and milk goats. Genomics. 2020;112(1):332–45.

    CAS  PubMed  Article  Google Scholar 

  16. Kulessa H, Turk G, Hogan BL. Inhibition of bmp signaling affects growth and differentiation in the anagen hair follicle. EMBO J. 2000;19(24):6664–74.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. Rishikaysh P, et al. Signaling involved in hair follicle morphogenesis and development. Int J Mol Sci. 2014;15(1):1647–70.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. Yue Y, et al. 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.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  19. Lin CM, et al. Long noncoding RNA expression in dermal papilla cells contributes to hairy gene regulation. Biochem Biophys Res Commun. 2014;453(3):508–14.

    CAS  PubMed  Article  Google Scholar 

  20. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018;172(3):393–407.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  21. Ma S, et al. Synchronous profiling and analysis of mRNAs and ncRNAs in the dermal papilla cells from cashmere goats. BMC Genomics. 2019;20(1):512–2.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  22. Zhao R, et al. Transcriptomic analysis reveals the involvement of lncRNA–miRNA–mRNA networks in hair follicle induction in Aohan fine wool sheep skin. Front Genet. 2020;11:590.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. Sulayman A, et al. Genome-wide identification and characterization of long non-coding RNAs expressed during sheep fetal and postnatal hair follicle development. Sci Rep. 2019;9(1):8501.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. Zhou G, et al. Integrative analysis reveals ncRNA-mediated molecular regulatory network driving secondary hair follicle regression in cashmere goats. BMC Genomics. 2018;19(1):222.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  25. Liu G, et al. Identification of microRNAs in wool follicles during anagen, catagen, and telogen phases in Tibetan sheep. PLoS One. 2013;8(10):e77801.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. Ren H, et al. Genome-wide analysis of long non-coding RNAs at early stage of skin pigmentation in goats (Capra hircus). BMC Genomics. 2016;17(1):67.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  27. Trapnell C, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. Pertea M, et al. Transcript-level expression analysis of RNA-seq experiments with HISAT. StringTie and Ballgown. 2016;11(9):1650.

    CAS  Google Scholar 

  29. Robinson, M.D., D.J. McCarthy, and G.K.J.B. Smyth, edgeR: a bioconductor package for differential expression analysis of digital gene expression data. 2010;26(1):139–140.

    Google Scholar 

  30. Benjamini Y, Yekutieli DJAOS. The control of the false discovery rate in multiple testing under dependency; 2001. p. 1165–88.

    Google Scholar 

  31. Karimizadeh E, et al. Analysis of gene expression profiles and protein-protein interaction networks in multiple tissues of systemic sclerosis. BMC Med Genet. 2019;12(1):199.

    CAS  Google Scholar 

  32. Wang J, et al. Identification of key pathways and genes related to the development of hair follicle cycle in cashmere goats. Genes (Basel). 2021;12(2):180.

    CAS  Article  Google Scholar 

  33. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. Kanehisa M, et al. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–d551.

    CAS  PubMed  Article  Google Scholar 

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

    CAS  PubMed  PubMed Central  Article  Google Scholar 

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

    CAS  PubMed  Article  Google Scholar 

  38. Statello L, et al. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol. 2021;22(2):96–118.

    CAS  PubMed  Article  Google Scholar 

  39. Shang F, et al. Expression profiling and functional analysis of circular RNAs in inner Mongolian cashmere goat hair follicles. Front Genet. 2021;12:678825.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. Gao Y, et al. KIF3C promotes proliferation, migration, and invasion of Glioma cells by activating the PI3K/AKT pathway and inducing EMT. Biomed Res Int. 2020;2020:6349312.

    PubMed  PubMed Central  Google Scholar 

  41. Bai WL, et al. LncRNAs in secondary hair follicle of cashmere goat: identification, expression, and their regulatory network in Wnt signaling pathway. Anim Biotechnol. 2018;29(3):199–211.

    CAS  PubMed  Article  Google Scholar 

  42. Wang W, et al. Long noncoding RNA: regulatory mechanisms and therapeutic potential in Sepsis. Front Cell Infect Microbiol. 2021;11(418):563126.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. Zhang X, et al. Mechanisms and functions of long non-coding RNAs at multiple regulatory levels. Int J Mol Sci. 2019;20(22):5573.

    CAS  PubMed Central  Article  Google Scholar 

  44. Sebastian-delaCruz M, et al. The role of lncRNAs in gene expression regulation through mRNA stabilization. Non-coding RNA. 2021;7(1):3.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  45. Mutz K-O, et al. Transcriptome analysis using next-generation sequencing. Curr Opin Biotechnol. 2013;24(1):22–30.

    CAS  PubMed  Article  Google Scholar 

  46. Yue Y, et al. Exploring differentially expressed genes and natural antisense transcripts in sheep (Ovis aries) skin with different wool Fiber diameters by digital gene expression profiling. PLoS One. 2015;10(6):e0129249.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. Stenn KS, Karnik P. Lipids to the top of hair biology. J Invest Dermatol. 2010;130(5):1205–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  48. Fan SM-Y, et al. Inducing hair follicle neogenesis with secreted proteins enriched in embryonic skin. Biomaterials. 2018;167:121–31.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. Jin M, et al. Molecular characterization and expression pattern of a novel keratin-associated protein 11.1 gene in the Liaoning cashmere goat (Capra hircus). Asian Australas J Anim Sci. 2017;30(3):328–37.

    CAS  PubMed  Article  Google Scholar 

  50. Khan I, et al. Mammalian keratin associated proteins (KRTAPs) subgenomes: disentangling hair diversity and adaptation to terrestrial and aquatic environments. BMC Genomics. 2014;15:779. https://doi.org/10.1186/1471-2164-15-779.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. Langbein L, et al. The catalog of human hair keratins. I. Expression of the nine type I members in the hair follicle. J Biol Chem. 1999;274(28):19874–84.

    CAS  PubMed  Article  Google Scholar 

  52. Gong H, et al. Identification of the ovine KAP11-1 gene (KRTAP11-1) and genetic variation in its coding sequence. Mol Biol Rep. 2011;38(8):5429–33.

    CAS  PubMed  Article  Google Scholar 

  53. Nachat R, et al. Peptidylarginine Deiminase isoforms are differentially expressed in the Anagen hair follicles and other human skin appendages. J Investig Dermatol. 2005;125(1):34–41.

    CAS  PubMed  Article  Google Scholar 

  54. Ahlawat S, et al. Skin transcriptome profiling of Changthangi goats highlights the relevance of genes involved in pashmina production. Sci Rep. 2020;10(1):6050.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. Brea-Calvo G, et al. Cell survival from chemotherapy depends on NF-κB transcriptional up-regulation of coenzyme Q biosynthesis. PLoS One. 2009;4(4):e5301.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  56. Lefebvre S, et al. Identification of Ectodysplasin target genes reveals the involvement of chemokines in hair development. J Investig Dermatol. 2012;132(4):1094–102.

    CAS  PubMed  Article  Google Scholar 

  57. Philpott MP, et al. Effects of interleukins, colony-stimulating factor and tumour necrosis factor on human hair follicle growth in vitro: a possible role for interleukin-1 and tumour necrosis factor-alpha in alopecia areata. Br J Dermatol. 1996;135(6):942–8.

    CAS  PubMed  Article  Google Scholar 

  58. Wu J, et al. Circannual rhythm in the skin gene expression of cashmere goat. bioRxiv. 2020:2020.04. 04.023044:23-44.

  59. Zouboulis CC, et al. Alterations in innate immunity and epithelial cell differentiation are the molecular pillars of hidradenitis suppurativa. J Eur Acad Dermatol Venereol. 2020;34(4):846–61.

    CAS  PubMed  Article  Google Scholar 

  60. Dai B, et al. The overexpression of tβ4 in the hair follicle tissue of Alpas cashmere goats increases cashmere yield and promotes hair follicle development. Animals. 2019;10(1):75.

    PubMed Central  Article  Google Scholar 

  61. Jara CP, et al. Glutamic acid promotes hair growth in mice. Sci Rep. 2021;11(1):15453.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  62. Yen C-LE, et al. A human skin multifunctional O-acyltransferase that catalyzes the synthesis of acylglycerols, waxes, and retinyl esters. J Lipid Res. 2005;46(11):2388–97.

    CAS  PubMed  Article  Google Scholar 

  63. Maier H, et al. Normal fur development and sebum production depends on fatty acid 2-hydroxylase expression in sebaceous glands. J Biol Chem. 2011;286(29):25922–34.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  64. Liu H, Wintour EM. Aquaporins in development–a review. Reprod Biol Endocrinol. 2005;3(1):1–10.

    Article  CAS  Google Scholar 

  65. Murphrey MB, Agarwal S, Zito PM. Anatomy, Hair. StatPearls. 2021. https://www.ncbi.nlm.nih.gov/books/NBK513312.

  66. Sasaki T, et al. Integrative analysis identifies activated anti-tumor immune microenvironment in lung metastasis of pancreatic cancer. Int J Clin Oncol. 2022;27(5):948–57.

    CAS  PubMed  Article  Google Scholar 

  67. Mullin J, et al. Transcription of canine toll-like receptor 2, β-defensin 1 and β-defensin 103 in infected atopic skin, non-infected atopic skin, healthy skin and the CPEK cell line. Vet Microbiol. 2013;162(2):700–6.

    CAS  PubMed  Article  Google Scholar 

  68. Mentlein L, et al. P014 FAM167A/BLK is a susceptibility locus in autoimmune diseases: characterisation of the FAM167 gene family. Ann Rheum Dis. 2018;77(Suppl 1):A19–9.

    Google Scholar 

  69. Wei JP, et al. EGFL6 expression in hair follicle central isthmus is dependent on the presence of terminal Schwann cells. Exp Dermatol. 2020;29(4):400–3.

    CAS  PubMed  Article  Google Scholar 

  70. Cheng C-C, et al. Hair follicle epidermal stem cells define a niche for tactile sensation. Elife. 2018;7:e38883.

    PubMed  PubMed Central  Article  Google Scholar 

  71. Katoh M. Function and cancer genomics of FAT family genes. Int J Oncol. 2012;41(6):1913–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  72. de Sousa N, et al. Hippo signaling controls cell cycle and restricts cell plasticity in planarians. PLoS Biol. 2018;16(1):e2002399.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  73. Klüppel M, et al. C4ST-1/CHST11-controlled chondroitin sulfation interferes with oncogenic HRAS signaling in Costello syndrome. Eur J Hum Genet. 2012;20(8):870–7.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  74. Kim YM, et al. SH3BP4 is a negative regulator of amino acid-rag GTPase-mTORC1 signaling. Mol Cell. 2012;46(6):833–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  75. Meng D, Frank AR, Jewell JL. mTOR signaling in stem and progenitor cells. Development. 2018;145(1):dev152595.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  76. Yanagisawa H, et al. Fibulin-5 is an elastin-binding protein essential for elastic fibre development in vivo. Nature. 2002;415:168.

    PubMed  Article  Google Scholar 

  77. Lee J, et al. Runx1 and p21 synergistically limit the extent of hair follicle stem cell quiescence in vivo. Proc Natl Acad Sci U S A. 2013;110(12):4634–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  78. Gao JG, Simon M. A comparative study of human GS2, its paralogues, and its rat orthologue. Biochem Biophys Res Commun. 2007;360(2):501–6.

    CAS  PubMed  Article  Google Scholar 

  79. Tashiro K, et al. GAREM, a novel adaptor protein for growth factor receptor-bound protein 2, contributes to cellular transformation through the activation of extracellular signal-regulated kinase signaling. J Biol Chem. 2009;284(30):20206–14.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  80. Kennedy SA, et al. Extensive rewiring of the EGFR network in colorectal cancer cells expressing transforming levels of KRAS(G13D). Nat Commun. 2020;11(1):499.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  81. Yamauchi K, Kurosaka A. Expression and function of glycogen synthase kinase-3 in human hair follicles. Arch Dermatol Res. 2010;302(4):263–70.

    CAS  PubMed  Article  Google Scholar 

  82. Frame S, Cohen P. GSK3 takes Centre stage more than 20 years after its discovery. Biochem J. 2001;359(1):1–16.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  83. Liu G, et al. Expression profiling reveals genes involved in the regulation of wool follicle bulb regression and regeneration in sheep. Int J Mol Sci. 2015;16(5):9152–66.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  84. Pangon L, et al. MCC inhibits beta-catenin transcriptional activity by sequestering DBC1 in the cytoplasm. Int J Cancer. 2015;136(1):55–64.

    CAS  PubMed  Article  Google Scholar 

  85. Rognoni E, Walko G. The roles of YAP/TAZ and the hippo pathway in healthy and diseased skin. Cells. 2019;8(5):411.

    CAS  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This research was funded by Xinjiang Autonomous Region Innovation Environment Construction Special Project (2021D04008, 2020Q035), Agriculture Research System of China (CARS-39-03), Innovation Project of Shandong Academy of Agricultural Sciences (13200214443101).

Author information

Affiliations

Authors

Contributions

CLW: Methodology, Writing-Original Draft, CKQ: Formal analysis, Investigation, XFF: Validation, Data Curation, XXH: Writing-Review & Editing, Conceptualization, KCT: Supervision, Project administration. All authors have read and approved the final version of the manuscript.

Corresponding authors

Correspondence to Xixia Huang or Kechuan Tian.

Ethics declarations

Ethics approval and consent to participate

Sample collection was carried out under license in accordance with the Guidelines for Care and Use of Laboratory Animals of China and all protocols were approved by the Animal Care and Use Committee of Xinjiang Academy of Animal Science (Approval number 2020008). All animal experiments were performed in accordance with the ARRIVE guidelines (https://arriveguidelines.org).

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

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

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wu, C., Qin, C., Fu, X. et al. Integrated analysis of lncRNAs and mRNAs by RNA-Seq in secondary hair follicle development and cycling (anagen, catagen and telogen) of Jiangnan cashmere goat (Capra hircus). BMC Vet Res 18, 167 (2022). https://doi.org/10.1186/s12917-022-03253-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12917-022-03253-0

Keywords

  • Goat hair follicles
  • Signaling pathways
  • Differentially expressed genes
  • Photoperiodism