The single-cell sequencing identifies 13 distinct cell populations in goat skins. To better define the cell types in Liaoning Cashmere Goat skins, we performed a single-cell RNA sequencing analysis in coarse type of Liaoning Cashmere Goat skins (CT-LCG) and fine type of Liaoning Cashmere Goat skins (FT-LCG). For each type goat, we obtained 551,655,085 and 556,307,136 reads. The sequencing saturation rates were 76.9% and 68.3%, respectively (Fig. 1a). We separated and sequenced 6,102 and 1,1197 cells, following strict quality control 3,443 and 9,608 cells were remained for downstream analysis in CT-LCG and FT-LCG.
We first defined different cell populations in cashmere goat skins, using the analysis of unsupervised clustering and t-distributed stochastic neighbor embedding dimension reduction (t-SNE). After standard quality control excluded cells including genes check <500 and mitochondrial genes coverage >10%, samples were detailed clustering analyzed by Seurat workflow definition 13 distinct population (Fig. 1b). The cells interaction between 13 clusters is shown in Fig. 1c. The average expression profile of the two samples is highly correlated, and the number of cells in each type is shown in Fig. 1d. We have identified the gene expression profile for per cluster, differential expression genes (DEGs) refer to the difference between the average expression of cells in clusters and that not in clusters (Fig. 1f). The top10 up expression genes of each cluster are shown in Table 1. In total, 2 kinds dermal papilla cells, 4 stem cells, 3 sheath cells, 1 keratinocyte, 1 fibroblast, 1 epidermal cell and 1 stromal cell identity was assigned in goat skin based on known cell type specific marker expression such as CD34, KRT19. We detected that genes enriched in distinct clusters differently marked various cell types. For instance, KRT14 and KRT5 are known expression markers in Epidermis stem cells (Cluster 1) (Fig. 1e). EDAR and DKK1 are well known to be highly expressed in Epidermis cells (Cluster 2). Dermal sheath cells (Cluster 6) were identified with specific expression of known marker gene BMP2. Genes enriched in the Cluster 7 strongly labeled in Dermal stem cells and expressed cell markers such as SOX2 and FOXM1. Hair follicle stem cells (Cluster 8, HFSCs) promote hair development, and were identified by multiple specific markers (e.g., KRT19; DSG3 and ITGA6). Skin stem cells (Cluster 12) also expressed selective genes, including CD34 and ITGB1. Finally, we found that cells from C10 showed characteristic of Dermal papilla cells (DPCs), which is known to be high expression level of CTNNB1, and C11 showed characteristic of Secondary hair follicle dermal papilla cells (SDPCs) such as high expression levels of VCAN and ACTA2. Moreover, Dermal fibroblasts cells markers (e.g., EGF and VEGFA), Inner root sheath cells markers (e.g., CDIP1; CDS1 and FOXN1), Outer root sheath cells markers (e.g., DOCK6 and KRT19), Stromal cells marker (e.g., ITGA6) and Keratinocytes cells marker (e.g., KRT6A) showed specificity expression in Cluster 0, Cluster 3, Cluster 5, Cluster 4 and Cluster 9, respectively. In short, these results emerged the comprehensive transcriptional heterogeneity of cell types, which make up cashmere goat skins and revealed molecular generate markers to make further efforts study in diversities and functions.
Importantly, the high expression genes of per cell group are shown in (Fig. 1g). KRT5 and KRT14 are highest expressed in Epidermis stem cell, and significantly enriched in Cluster 0, Cluster 2, Cluster 7 and Cluster 8. CCL27, ADIRF and LY6D had obvious difference expression in Epidermis stem cells. KRT10 are highest expressed in Cluster 4 and Cluster 5, LOC102179515, NFIB, BGN and SFPRP1 was overexpressed on Cluster 8. Interestingly, CXCL14 was overexpressed on Cluster 8 and Cluster 11, which may be related to hair development. Indeed, KRTDAP and CST6 are highest expressed in Cluster 9. Notably, KRT35 and KRTAP11-1, may be related to cashmere character12, whereas in goat their expressions are confined to Dermal papilla cells. Additionally, LEF1 is considered to be a hair-inducing gene of dermal papilla cells13. While the expression of the positive genes, IGFBP5 and IGFBP6, are discovered highly expressed in SDPCs. ALX4, which is known may be related to cashmere fineness14, only expressed in DPCs and SDPCs. CCL21, GNG11 and TFPI2 are highest expressed in HFSCs, VIM and CLEC3B were overexpressed on DPCs and HFSCs. KRT19 and SOX9 were proven to exist in the anagen secondary HFSCs15, which have the high expressed in our results in HFSCs. Some genes are only expressed in certain cells, such as PLVAP, PECAM1 and CENPF (Fig. 1h). These significant up expressed genes, especially over-expressed in each functional cluster, and may be functional candidate genes associated with fibrogenesis. Collectively, these data established the classification of skin cells and the selective genes expression pattern of cashmere goat.
Table 1. The top 10 genes in per Clusters.
Expression
|
Cluster 0
|
Cluster 1
|
Cluster 2
|
Cluster 3
|
Cluster 4
|
Cluster 5
|
Cluster 6
|
Cluster 7
|
Cluster 8
|
Cluster 9
|
Cluster 10
|
Cluster 11
|
Cluster 12
|
Top1
|
KRT10
|
KRT14
|
KRT14
|
LOC108634769
|
KRT10
|
KRT10
|
KRT1
|
KRT14
|
KRT5
|
CST6
|
RPS12
|
LOC108634769
|
TMSB4X
|
148.51
|
334.64
|
170.19
|
232.43
|
232.68
|
322.73
|
197.52
|
194.81
|
125.12
|
792.65
|
132.77
|
82.62
|
102.97
|
Top2
|
KRT14
|
KRT5
|
KRT5
|
KRT10
|
LOC102177275
|
KRT1
|
KRT10
|
KRT5
|
KRT14
|
KRTDAP
|
LOC102168701
|
LOC102168701
|
LOC108634769
|
127.88
|
158.27
|
129.94
|
173.85
|
179.90
|
197.82
|
155.68
|
112.26
|
112.85
|
300.61
|
118.94
|
73.91
|
97.50
|
Top3
|
LOC102177275
|
LGALS7
|
LOC102177275
|
KRT1
|
KRT1
|
S100A7A
|
LOC108634769
|
PTMA
|
LOC102179515
|
KRT10
|
RPS8
|
COL1A1
|
VIM
|
108.49
|
108.19
|
110.33
|
128.36
|
132.21
|
155.10
|
116.21
|
95.76
|
106.09
|
265.48
|
115.23
|
71.50
|
97.21
|
Top4
|
CST6
|
LOC102168701
|
RPS8
|
RPLP1
|
LGALS7
|
KRTDAP
|
RPLP1
|
RPS8
|
LOC102168701
|
SBSN
|
RPLP1
|
RPS8
|
TPT1.1
|
101.48
|
105.19
|
94.81
|
93.46
|
110.19
|
99.43
|
71.73
|
90.80
|
93.90
|
174.83
|
97.81
|
70.63
|
84.12
|
Top5
|
KRT5
|
TPT1.1
|
RPLP1
|
LOC102168701
|
ACTG1
|
SPINK7
|
LOC108635081
|
LOC102168701
|
RPS8
|
DMKN
|
RPLP0
|
VIM
|
LOC102168701
|
88.78
|
101.40
|
80.37
|
90.56
|
103.55
|
97.56
|
70.77
|
87.07
|
91.65
|
109.14
|
87.65
|
62.30
|
80.81
|
Top6
|
ACTG1
|
RPLP1
|
LOC102168701
|
RPS8
|
CST6
|
LGALS7
|
LGALS7
|
RPLP1
|
RPS12
|
KRT1
|
RPS2
|
RPLP1
|
RPS8
|
81.70
|
98.20
|
79.72
|
87.21
|
88.81
|
86.67
|
66.77
|
85.31
|
83.15
|
107.72
|
80.94
|
61.71
|
78.59
|
Top7
|
KRT1
|
RPS8
|
TPT1.1
|
RPS12
|
KRTDAP
|
RPLP1
|
JUN
|
RPS12
|
TPT1.1
|
HSPB1
|
RPL32
|
S100A4
|
LOC102182395
|
70.62
|
91.93
|
70.67
|
80.54
|
84.74
|
83.42
|
61.55
|
67.79
|
79.14
|
89.67
|
76.40
|
61.11
|
75.92
|
Top8
|
LGALS7
|
S100A2
|
RPS12
|
S100A7A
|
RPS8
|
LOC102168701
|
KRTDAP
|
LGALS7
|
RPLP1
|
LGALS7
|
LOC102184223
|
RPS12
|
RPS12
|
70.56
|
90.87
|
68.12
|
67.47
|
80.17
|
79.57
|
58.47
|
66.43
|
76.50
|
81.33
|
75.60
|
60.18
|
70.72
|
Top9
|
LOC102168701
|
RPS12
|
LGALS7
|
LGALS7
|
HSPB1
|
SBSN
|
RPS8
|
RPS19
|
LOC102177275
|
LOC102177275
|
RPS5
|
TPT1.1
|
LOC108635081
|
67.65
|
84.82
|
66.98
|
66.95
|
75.44
|
76.47
|
55.76
|
64.08
|
73.90
|
80.80
|
74.67
|
57.61
|
67.08
|
Top10
|
RPLP1
|
RPS19
|
RPS19
|
LOC102169125
|
RPLP1
|
RPS12
|
UBB
|
RPS2
|
RPS19
|
RPLP1
|
RPL31
|
RPS24
|
RPLP1
|
66.06
|
80.19
|
65.51
|
64.82
|
74.39
|
70.96
|
55.69
|
60.15
|
72.41
|
70.05
|
73.46
|
50.28
|
66.72
|
Pseudo-Timeline analysis revealscelldevelopment in cashmere goat skins. To further investigate the classification of skin cell populations and developmental relationship, we perform cell trajectory analysis using the Monocle analysis toolkit16. To enable a single-cell to differentially expressed genes during development based on their trajectories. Cells were calculated by Monocle 2 with an unsupervised way via maximizing transcription similarity between consecutive cell pairs. Therefore, this method can be used for identify various phases of skin development in cashmere goats. Monocyte 2 has a large overlap in each cluster, and the subpopulation of cell classification is distributed in the whole skin, indicating versatility.
To explore which gene regulates the cashmere fiber development progression through per axis, we classified and clustered the gene expression changes with pseudo-time. Pseudo-time results reveal the skin development of cashmere goat is divided into 15 state, every kind of cell is involved in almost all the development process (Fig. 2a, 2b). In support of the draw order in pseudo-time trajectory, our analysis revealed the organizes of genes that performed differentially expressed in each branch, which closely related to the genes that are famous in participated in cashmere growth and hair development. Interestingly, KRTAP11-1, TCHH, KRTDAP and S100A8 that were briefly up-regulated and then down-regulated were observed. However, a subset of genes such as KRT5, RPS7, RPL3 and RPL5 were sustained up-regulated over time (Fig. 2c). Further confirming the validity of pseudo-timeline analysis, the top 6 genes such as KRTAP11-1 and TCHH were enriched (Fig. 2d), which GO analysis enriched Intermediate filament and Keratinization are important pathways for cashmere growth, indicating these genes may relate to cashmere formation and development. KEGG pathways of these genes mainly enriched to WNT, TGF beta and PI3K-Akt signaling pathway. These data and analyses revealed the clear developmental stages of goat skins, the differential genes expressions along the time line of bifurcations indicates that different molecular paths guide the development of different populations.
Identification and analysis DPCs andSDPCs in cashmere goat skins. Previous work showed that Dermal papilla cells can induce hair follicle regeneration, determine the size of hair follicle to a certain extent, and play an essential role in hair follicle cycle change process. Dermal papilla cells subtly grasp the continuous growth of mammalian hair, including cashmere. There are 100 high expression genes in DPCs, and the top10 genes exhibited in Fig. 3. LOC102184693, KRT35, KRTAP11-1, KRTAP3-1 and MT4 almost only expressed in DPCs. In hair follicle and goat skin, KRT35 and KRTAP11-1 are frequently appeared. Moreover, we detected two crucial activators HOXC8 and RSPO1 were enriched in DPCs, as reported before17. In addition, the family of RPS and RPL are significantly expressed in DPCs, such as RPS8, RPS2, RPL0 and RPL32. RPS12, as the key effector of cell competition caused by other RP gene mutations, has the highest expression in DPCs. Among them, PTGER4 and ESR1 are more abundant in DPCs consistent with previous studies18, 19. To a better understanding the functioning approaches of differentially expressed genes in DPCs, GO and KEGG enrichment analysis are performed. Finally, 452 pathways were enriched such as keratin filament and hair follicle morphogenesis. By KEGG analysis some pathways related to cashmere growth and development were enriched, including TGF-beta signaling pathway, Alcoholism, PPAR and Notch signaling pathway.
Recently, research in this field has indicated that hair type can be determined by SDPCs. In SDPCs, the top10 up-expression genes exhibited in Fig. 4a, 4b. Of note, COL1A1, COL1A2, COL3A1, S100A4, VIM, CRABP2, CCDC80, ACTG2, CCL26 and CXCL12 were solely abundant in SDPCs. COL1A2, COL3A1 and CCDC80 were reported to be differentially expressed in different cashmere fibers12. Furthermore, the expression of S100A4 was the highest in SDPCs, which is may involve in Yangtze River Delta white goat hair growth 20. Also, NFKBIA had a highly enriched in SDPCs, which may play an important role in cashmere growth12, 21. Similar, CRABP1, a constant marker of DPCs, which has a high expressed in SDPCs, and expresses throughout the whole stages of hair cycling22. In order to further excavate the critical genes regulating cashmere fineness, enrichment was analyzed in SDPCs. A total of 848 GO terms were significantly enriched, as expected, hair follicle morphogenesis, hair follicle development and intermediate filament cytoskeleton were significantly-enriched GO terms. The top 20 differentially expressed GO terms were showed in Fig. 4c. The top 20 pathways by KEGG analysis were significantly enriched, such as PI3K-Akt, TGF beta and WNT signaling pathway (Fig. 4d). BMP3, BMP4, BMP7, FGF7, FGF2, COL1A1, IGF1, and WNT5A enriched these pathways.
To explore difference genes in development period, we took a closer analyzed of Pseudo-Timeline Analysis of DPCs and SDPCs, respectively. Trajectory analysis provided the development process and five differentiation state in DPCs (Fig 5a, 5b). CST6, KRT10, KRTDAP and DMKN decreased first, then rise and then fall at DPCs (Fig 5c, 5d), and yet KRT35 and MT4 were momentary upregulated and then downregulated. SOX4 and CEBPB continually upregulated and KRT14 and KRT1 downregulated. Pseudo-time analysis demonstrated a total of 5 state in SDPCs development stage (Fig 6a, 6b). CXCL8 and CD34, after a transiently period of stability, it began to downregulated (Fig 6c). Notably, CXCL4, COL1A1, COL1A2 and COL3A1 were continually downregulated and SRGN, NFKB1 and RGCC showed late specific expression (Fig 7d). As expected, GO analysis showed that DPCs and SDPCs both significantly-enriched KRT14,CXCL4 and COL1A1 in the development process, and with the development of pseudo-timeline on the axis, these genes involved in the development were down regulated. KRT35, SRGN and NFKB1 showed specific expressions in development later.CXCL4, COL1A1, KRT14 and KRT35 may be the key candidate genes for regulation DPCs and SDPCs.
Weighted Gene Correlation Network Analysis (WGCNA) of cashmere goat skins. In order to further identify the genes related to cashmere growth, we conducted the WGCNA on cashmere goat skins. WGCNA, a relatively new statistical analysis, is usually used to construct networks based on gene correlations and recognize intramodular hub genes23. The cluster dendrogram consists of 31 co-expression modules, including tan, greenyellow, white, green and other colors (Fig 7a). We also found that compared with others, the blue module has a strong co-expression relationship by building a topological overlap matrix (TOM) diagram (Fig. 7b). The darker the color, the higher the topological overlap.
As shown in Fig. 7c, all genes in 13 type cells were selected to produce the module trait relationship. Interestingly, before there is almost no similarity between each module, suggesting that most modules had certain promoting effect on each cluster. The significantly enrichment modules in secondary hair follicles dermal papilla cells is blue, in dermal papilla cells is green, and in stromal cell is gray60, and the number of genes in per module was 2629, 1839, and 256, respectively. The blue module of SDPCs cluster related to cashmere growth has the deepest color, suggesting that genes in this module may have a positive effect on cashmere fineness. We identified hub gene COL1A1, COL1A2, CCL26, CCDC80, HOXA13, CXCL8, KRT24, and BMP3 has significant connectivity, which may associate with cashmere development.
Enrichment analysis was used to explore the biological functions of genes in these modules, to further determine the relationship between these genes and cashmere growth. The GO analysis results of the core genes in blue module were displayed in Fig. 7d. It was found that these genes are suggested to cytoplasm, integral component of the membrane.
The detailed network information of core genes for 4 selected modules including blue, brown, green and gray60 modules was showed in Fig 8. In the gene networks, some hub genes interact with many other genes, which indicates that they are more likely to be necessary. COL1A1, CCDC80 and HOXA13 are associated with many genes in blue module. CXCL8, NFKBIE and CCL24 are actively in green module. KRT3, KRT4 and KLK12 are positively linked to other genes in gray60 module. KRT19, KRT38, NFIB and LHX2 are positive genes in brown module. COL1A1, with a significant connectivity in SDPCs, may serve as a target for cashmere formation. IF staining with COL1A1 showed that COL1A1 is expressed in SDPCs, which is as predicted by scRNA-seq and WGCNA (Fig. 9).
Analysis and verification of genes difference in SDPCs of CT-LCG and FT-LCG. A total of 16601 genes were found in SDPCs of two samples (Fig. 10a), there are 15 representative DEGs between CT-LCG and FT-LCG exhibited in Fig. 10b. TCHH, S100A8 and IGF2 were highly expressed in CT-LCG, nevertheless, CD74, CCL27, COL1A1 and CXCL8 were highly expressed in FT-LCG. In order to further excavate the critical genes regulating cashmere fineness, enrichment was analyzed in SDPCs. There are 48 and 45 significant pathways in CT-LCG and FT-LCG. CT-LCG are mainly enriched to Jak-STAT signaling pathway, TGF-beta signaling pathway and PI3K-Akt signaling pathway (Fig. 10c), and FT-LCG also enriched in TGF-beta signaling pathway (Fig. 10d). COL1A1 enriched to PI3K-Akt signaling pathway and CXCL8 enriched to NF-kappa B signaling pathway, which had been reported may be related to cashmere.
We identified ACTA2 (α-SMA) as a marker gene of SDP cells by immunofluorescence, α-SMA protein was almost exclusively expressed on secondary hair follicles, as shown in Fig. 11. To further determine the relationship between these genes and cashmere fineness, we validated COL1A1 and CXCL8 in SDPCs by immunofluorescence assay. We have made immunofluorescence test on the skin of three coarse and three fine Liaoning Cashmere Goats. The differential expression of COL1A1 in coarse and fine type skin indicates that the number of SHFs in fine skin is more distributed, arranged tightly and with smaller cashmere diameter, while the number of SHFs in coarse skin is less distributed, arranged sparsely and with larger cashmere diameter (Fig. 12). Wang have shown that COL1A1 has difference in CT-LCG and FT-LCG24. In the results of immunofluorescence, the nucleus is blue and the target protein is green. PHF is the primary hair follicle, SHF is the secondary hair follicle, which is the place where cashmere grows. SeG is sebaceous glands. They form the hair follicle cluster. The SDPCs exist in the secondary hair follicles. CXCL8 is the differential expression gene of CT-LCG and FT-LCG screened by us. The result shows that the expression level of CXCL8 (IL-8) protein in CT-LCG is higher than that in FT-LCG (Fig. 13). From the results, we can see that the results of immunofluorescence are consistent with the results of single-cell sequencing.
The OD value of protein standard was measured at the wavelength of 568 nm, and the standard curve was established. The regression equation and variance coefficient were y = 0.6965x + 0.0141 and R2 = 0.997. The results showed that the expression level of ACTA2 (α-SMA) and CXCL8 in FT-LCG was lower than that in CT-LCG, as shown in Fig. 14.