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.

Keywords: Differentially expressed genes; Goat hair follicles; Photoperiodism; Signaling pathways.

Fig. 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

Fig. 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 (log₂ (fold change) ≥1 and p-value< 0.05)

Fig. 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

Fig. 4

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

Fig. 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

Fig. 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

Fig. 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

Fig. 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 (log₂ (fold change) ≥1 and p-value< 0.05)

Fig. 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

Fig. 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

Fig. 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

Fig. 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

Fig. 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

Fig. 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

Fig. 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

Fig. 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

Fig. 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)