Next Article in Journal
A Breeding Plumage in the Making: The Unique Process of Plumage Coloration in the Crested Ibis in Terms of Chemical Composition and Sex Hormones
Previous Article in Journal
Inhibition of Cell Apoptosis by Apicomplexan Protozoa–Host Interaction in the Early Stage of Infection
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Animals
Volume 13
Issue 24
10.3390/ani13243818
animals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Interpretation of the Yak Skin Single-Cell Transcriptome Landscape

by
Qingbo Zheng
1,2,†,
Na Ye
1,2,†,
Pengjia Bao
1,2,
Tong Wang
1,2,
Chaofan Ma
1,2,
Min Chu
1,2,
Xiaoyun Wu
1,2,
Siyuan Kong
3,
Xian Guo
1,2,
Chunnian Liang
1,2,
Heping Pan
4,* and
Ping Yan
1,2,5,*
1
Key Laboratory of Animal Genetics and Breeding on Tibetan Plateau, Ministry of Agriculture and Rural Affairs, Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences, Lanzhou 730050, China
2
Key Laboratory of Yak Breeding Engineering of Gansu Province, Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences, Lanzhou 730050, China
3
Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen 518120, China
4
Life Science and Engineering College, Northwest Minzu University, Lanzhou 730030, China
5
Institute of Western Agriculture, The Chinese Academy of Agricultural Sciences, Changji 831100, China
*
Authors to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Animals 2023, 13(24), 3818; https://doi.org/10.3390/ani13243818
Submission received: 2 November 2023 / Revised: 5 December 2023 / Accepted: 7 December 2023 / Published: 11 December 2023
(This article belongs to the Section Animal Genetics and Genomics)
Browse Figures
Versions Notes

Abstract

:

Simple Summary

Yaks can grow in extremely harsh natural environmental conditions, such as high altitude, low temperature, and hypoxia, and hair plays an important role. As an appendage of the skin, hair follicles involve the specialization of various types of cells and intercellular communication during morphogenesis. Here, we used single-cell RNA sequencing technology to identify 11 cell types from the scapular skin of yaks, constructed a transcription map of the main cells of hair follicles in the skin, and described the heterogeneity of DP cells and dermal fibroblasts in the hair follicle cycle. Our findings provide a molecular landscape of the fate determination process of dermal cell lineages. Our research provides valuable resources for further exploration of the molecular pathways involved in hair follicles.

Abstract

The morphogenesis of hair follicle structure is accompanied by the differentiation of skin tissue. Mammalian coats are produced by hair follicles. The formation of hair follicles requires signal transmission between the epidermis and dermis. However, knowledge of the transcriptional regulatory mechanism is still lacking. We used single-cell RNA sequencing to obtain 26,573 single cells from the scapular skin of yaks at hair follicle telogen and anagen stages. With the help of known reference marker genes, 11 main cell types were identified. In addition, we further analyzed the DP cell and dermal fibroblast lineages, drew a single-cell map of the DP cell and dermal fibroblast lineages, and elaborated the key genes, signals, and functions involved in cell fate decision making. The results of this study provide a very valuable resource for the analysis of the heterogeneity of DP cells and dermal fibroblasts in the skin and provide a powerful theoretical reference for further exploring the diversity of hair follicle cell types and hair follicle morphogenesis.

1. Introduction

Hair is a highly keratinized tissue produced by mammalian hair follicles [1]. Hair follicles are a kind of skin appendage with independent functions and periodic growth and can be divided into epithelial cell groups and dermal cell groups according to cell origin and function [2,3,4]. The hair follicle is composed of a connective tissue sheath (CTS), hair bulb, inner root sheath (IRS), outer root sheath (ORS), and hair shaft (HS), and the HS is located in the center of the hair follicle and consists of three layers of cells: hair cuticle, medulla, and cortex [3,4,5]. The interaction between the epithelial cell group and dermal cell group contributes to the initial stage of hair follicle development [6]. Hair follicles undergo a continuous cycle of growth, that is, the cycle of telogen, anagen, and catagen, and depend on different cell signal exchanges between the epidermis and dermis [7,8]. Various types of cells in the hair follicle structure can control the proliferation and differentiation of hair follicle epidermal cells by secreting many signaling molecules, such as growth factors and cytokines.
Dermal fibroblasts are the most abundant cell type in the dermal layer of the skin, and their main function is to secrete the components of the extracellular matrix (ECM) to provide structural support [9]. Fibroblasts in the dermis are divided into different morphological and functional subgroups according to their spatial localization [9,10,11]. Fibroblasts connected to the epidermal basement membrane in the dermis are called papillary fibroblasts (PFs), and those located in the lower reticular dermis fibroblasts are called reticular fibroblasts (RFs). Driskell et al. [11] used lineage tracing to trace the origin of fibroblasts at different levels of the dermis in a mouse embryonic model and found that dermal fibroblast progenitor cells differentiated into the PF lineage and RF lineage at 16.5 weeks in mouse embryos. The process of skin aging mainly affects the dermis, the spinous process structure at the junction of the epidermis and dermis is significantly reduced, and the dermis shrinks [12,13,14]. The lineage relationship between all fibroblast subtypes in the dermis plays an important role in development, homeostasis, aging, wound repair processes, and maintaining the healthy state of the dermal microenvironment [11,15,16]. In recent years, single-cell RNA sequencing (scRNA-seq) has provided unprecedented insights into the specific transcription profiles of cell subsets in the skin, particularly regarding fibroblasts and their relationship with other cells in the dermis. Guerrero-Juarez et al. [17] used scRNA-seq technology to explore the diversity of skin wound fibroblasts and obtained 12 fibroblast clusters, some of which may be in a continuous differentiation state towards the contractile phenotype; some cell clusters may represent different fibroblast lines, and some subsets of cell lines express hematopoietic marker genes. Vorstandlechner et al. [18] used scRNA-seq to analyze the heterogeneity of human skin fibroblasts and identified six FB subgroups and found that the newly identified FB subgroups did not overlap with the markers commonly used to identify papillary and reticular fibroblasts. Thompson et al. [19] studied the accessibility characteristics of myofibroblast markers in neonatal fibroblasts and showed that different fibroblast subtypes, such as PF and RF, can support myofibroblast status.
Dermal papilla (DP) cells control the formation and size of hair follicles and are the control center of hair follicle growth and the hair follicle cycle. The DP is considered the central player in hair follicular cycle regulation, being triggered by paracrine signaling, and it serves as a reservoir for multipotent stem cells and various growth factors, participating in the regulation of hair follicle development and growth [20]. The number of DP cells in hair follicles and the size of the dermal papilla area formed by them are key to distinguishing different types of hair follicles [21]. The number of DP cells in different hair follicle cycles affects the occurrence of hair follicle morphology [22]. During the late anagen phase, adjacent cells to DP cells undergo apoptosis, causing the hair follicle cycle to enter catagen; the hair follicle enters telogen, while DP cells remain in dense spherical form; as the follicle cycle transitions back to anagen, DP cells become loose, engulfed by new hair bulbs and become slender [5]. Rompolas et al. [23] removed DP cells from mouse skin hair follicles during telogen via two-photon confocal microscopy and found that hair follicles that removed DP cells could not enter anagen, while hair follicles with intact DP cells could enter anagen normally. When the number of DP cells decreased, the process of entering anagen in hair follicles was delayed.
The yak is mainly distributed in the Qinghai-Tibet Plateau and its surrounding areas. It is a valuable livestock and poultry variety resource in the Qinghai-Tibet Plateau [24]. Yaks have thick skin and dense hair, which is one of the important reasons for their ability to adapt to cold and strong-ultraviolet environments. The yak has a unique mixed hair structure; its hair is composed of coarse hair, two-type hair, and fluff, and yak fluff is also a precious textile material [25]. The hair follicles of the yak show obvious seasonal periodic changes, accompanied by the growth and shedding of villi. Hair follicle morphogenesis is the result of the interaction of many cell types. Compared with other common villus-producing animals, there are few studies on different cells during the hair follicle cycle development of the yak. The study of the heterogeneity of dermal cells in the telogen and anagen phases of hair follicles will help to reveal the mechanism of dermal cell interaction and signal communication, which is significant for understanding the underlying processes in hair follicle biology.
To systematically understand the diversity of skin dermal cells, and to fully understand the pedigree heterogeneity of DP cell and dermal fibroblasts in the dermis and the fate regulatory mechanism of different cell subtypes, here, we used scRNA-seq technology to identify 18 cell clusters and 11 cell types from the scapular skin of the yak, constructed the main cell transcription map in hair follicles, and focused on the regulatory mechanism of fate determination in the specialization process of different fibroblast subsets. This work not only enriches people’s understanding of the morphogenesis of hair follicles in different cycles of yak, but also clarifies the differences in fibroblast subsets and their mechanism of action in the maintenance of skin homeostasis and provides a theoretical reference for further exploring the molecular mechanism of hair follicle cycle development in the yak.

2. Materials and Methods

2.1. Experimental Animals

The yak skin samples used in this study were all derived from the Tianzhu White Yak Breeding Farm in Wuwei City, Gansu Province. We selected three representative female yaks with a body weight of 2–3 years old consistent with the average body weight of their breeds and good health. According to the characteristics of yak hair follicle cycle development, we collected skin tissue of the yak scapula for the preparation of single-cell suspension during hair follicle anagen and telogen [1]. Before collection, the hair on the shoulder blade of the yak was trimmed by using elbow scissors, 2% lidocaine was injected subcutaneously for local anesthesia, and then the skin tissue was transferred using a skin sampler. From the collected skin, we removed as much subcutaneous fat tissue as possible, disinfected the sample with 75% alcohol, and washed with D-PBS to remove the alcohol. After treatment, the tissue was immediately stored in tissue protection solution and brought back to the laboratory for suspension preparation, and the other part was placed in paraformaldehyde for verification experiments.

2.2. Single-Cell RNA Sequencing of Yak Hair Follicles

The preparation of single-cell suspension and scRNA-seq analysis were performed according to the previous report [26]. After mixing the cell suspension of three samples into a mixture, it can be used for sequencing when the cell survival rate is ≥85%. The formation of single-cell GEM was performed using 10 × Genomics’ Chromium Single Cell 3’v3 Gel Beads kit. According to the instructions of the Chromium Single Cell 3’v3 Gel Beads kit, PE150 sequencing was performed using the Illumina NovaSeq 6000 platform. STAR in CellRanger was used to compare the readings with the splicing recognition of the reference genome. After CellRanger was completed, the next step is to cluster the scRNA-seq data. In this study, the “filtered_gene_bc_matrices” file generated via CellRanger analysis was read using Seurat in R environment (R version: 4.0.5), and then low-quality cells were filtered using the FilterCells function for different samples. The NormalizeData function was used to standardize the data, and then the expression matrix was subjected to dimensionality reduction clustering and cell type identification.

2.3. Construction of the Cell Differentiation Trajectory

To further analyze specific cell clusters and construct the differentiation trajectory of DP cells and fibroblast lines, we first used Seurat’s “SubsetData” function to extract fibroblast cell lines and DP cell clusters. Monocle (v2.18.0) was used to sort individual cells along pseudotime according to the official tutorial. To perform quasitime sorting on specific cell types, we use the “newCellDataSet” function in Monocle to construct Monocle objects. The state was set according to the cell cluster label identified using seurat, and the branch-specific expression genes were prepared using the “BEAM” function. The heat map was drawn using Monocle’s “plot_genes_branched_heatmap” function, and the genes with qval < 10−4 were regarded as input genes, after which the gene clusters were further divided into four clusters according to the k-means.

2.4. GO Enrichment Analysis

With the help of Metascape (http://metascape.org/gp/index.html#/main/step1) (accessed on 26 June 2023) online software, GO analysis of characteristic genes was performed. The gene names were transformed into gene IDs through the R package org.HS.eg.db (v3.12.0) annotation package, and the possible pathways were enriched with the clusterProfiler (v3.18.1), ggplot2 (v3.3.5), and enrichplot (v1.10.2) packages.

2.5. Cell Communication Analysis

CellChat (v1.5.0) was used to analyze intercellular communication. CellChat calculates the communication probability at the signaling pathway level by summarizing the communication probability of all ligand–receptor interactions associated with each signaling pathway. A circle plot was used to show the number and strength of interactions between cells. A bubble diagram was used to display the interaction ligand–receptor pairs, and Seurat was used to plot the gene expression distribution of signal genes related to the signaling pathway.

2.6. Immunofluorescence Analysis

The skin tissue was taken out from the 4% paraformaldehyde fixative. The skin tissue was trimmed with a scalpel, and the dehydration box was placed in a dehydrator (Diapath, Martinengo, BG, Italy) for dehydration and waxing leaching with different concentrations of gradient alcohol (Sinopharm Chemical ReagentCo., Ltd., Shanghai, China) and xylene (Sinopharm Chemical ReagentCo., Ltd., Shanghai, China). Then the melted wax was placed in the embedding box, and the tissue was taken out from the dehydration box and placed in the embedding box before the wax solidified and cooled at −20 °C. The tissue block was cut into 4 μm slices using a slicing machine (Shanghai Leica Instrument Co., Ltd., Shanghai, China), and then the slices were unfolded in 40 °C warm water. The tissue was fished up with a glass slide and baked at 60 °C. The sections were dewaxed in xylene, alcohol, and distilled water to water, and EDTA antigen retrieval buffer (Servicebio, Wuhan, China) was used for antigen repair. After natural cooling, the sections were placed in PBS (Servicebio, Wuhan, China) for washing. During the immunofluorescence single-labeling staining experiment, the slides were incubated overnight with the primary antibody at 4 °C, then placed in PBS and washed, and the corresponding second antibody was added and incubated in darkness at room temperature for 50 min. In immunofluorescence homologous double-labeling staining experiments, the first primary antibody needs to be added. After incubation at 4 °C overnight, the HRP-labeled secondary antibody of the corresponding species was added to cover the tissue and incubate at room temperature for 50 min. After washing the glass slides, TSA was added and incubated in the dark for 10 min, and then heated after washing. After that, the second primary antibody was added and incubated overnight at 4 °C. After washing, the fluorescent secondary antibody of the corresponding species of the primary antibody was added and incubated for 50 min. After incubation with primary and secondary antibodies, the slides were nuclear-stained with DAPI (Servicebio, Wuhan, China), autofluorescence quencher (Servicebio, Wuhan, China) was added for 5 min, and the sections were rinsed with running water for 10 min. Finally, the slides were examined using a fluorescence microscope (Nikon, Tokyo, Japan). The primary antibodies MSX1, VIM, and DLX3 were from Bioss (Bioss, Beijing, China).

3. Results

3.1. Single-Cell Data Set Quality Statistics

Using scRNA-seq, 26,573 single cells were obtained from the scapular skin of yaks. The data quality of each sample is shown in Table 1. The valid barcode of two samples is higher than 95%, and the genome reads mapped is greater than 85%. The number of cells detected in the telogen sample was 14,573, and the number of genes detected was 17,796; the number of cells detected in the anagen sample was 12,000, and the number of genes detected was 18,420. The overall data quality is at an ideal level.

3.2. Identification of Main Cell Types and Analysis of Characteristic Gene Expression

To construct the transcription map of yak skin tissue, we integrated the hair follicle telogen and anagen data, created a gene cell expression matrix, and reduced the dimension via UMAP clustering to determine 18 clusters (Figure 1a). We obtained 18 clusters of cell types. According to the comparison of the marker genes in the published papers with the specific genes expressed in each cluster, 11 cell types were identified. It was found that IFE-DC and DP cells were highly expressed during anagen (Figure 1b). Furthermore, cluster 1 and 11 (interfollicular epithelium differentiation cell, IFE-DC) marker genes KRT1 [27], KRT10 [27], and SBSN [27] were also found. Cluster 2 expressed the UHF cell marker genes KLK10 [28] and KRT79 [28]. Cluster 4 expressed the HS cell marker genes IGFBP5 [29] and LHX2 [30]. Cluster 5 expressed the DP cell marker genes LEF1 [31], MSX1 [32], and HOXC13 [32]. Cluster 7 expressed the epidermal cell lineage marker genes KRT5 [27] and KRT15 [33,34]. Cluster 8 expressed hair follicle stem cell (HFSC) cell marker genes SOX9 [35,36,37,38] and TCF4 [38]. Cluster 9 expressed the IRS cell marker genes KRT71 [39] and KRT25 [40]. Figure 1c shows the expression of some marker genes. Cluster 6 expressed dividing the fibroblast (DF) cell marker genes WNT6 [19] and NFIC [19]. Cluster 13 expressed the PF cell marker genes LRIG1 [11] and SCEL [19]. Cluster 14 expressed the RF cell marker genes LIPH [19], MGST1 [19], NEXN [13], and TPM1 [13]. Cluster 16 expressed the myofibroblast marker genes TGFBR2 [19] and FBLN2 [17]. Figure 1d shows the expression of some marker genes in different cells during anagen and telogen. The expression of FBLN2 in myofibroblasts was mainly in hair follicle anagen.
After identifying the cell types of yak hair follicles, we observed expression changes in some characteristic genes in hair follicles during telogen and anagen, and found that genes such as STAT3 and TCF4 were upregulated during anagen (Figure 2). Current studies have found that SHH, WNT/β-catenin/LEF1, and STAT3 signals are activated during the transition of hair follicles from telogen to anagen [3,41,42], indicating that STAT3 and other genes play a very important role in the transition of the hair follicle cycle in early anagen. The expression level of RSPO1 was significantly upregulated during the transition from telogen to anagen in mouse dorsal hair follicles [43]. We found that RSPO1 was expressed in yak hair follicles during hair follicle anagen and telogen. In a study of humans, it was found that RSPO1 was mainly expressed in dermal fibroblasts, and the expression of RSPO1 was not detected in cultured keratinocytes [44,45]. The expression level of TGFBR2 in myofibroblast cells was increased during the transition from telogen to anagen. Furthermore, we found that WNT3A was also up-regulated during hair follicle anagen. Guo et al. used qRT-PCR and Western blotting to detect the telogen and anagen of mouse hair follicles and found that WNT3A had the highest expression in the anagen of hair follicles and that the expression in telogen was weak [46].

3.3. Fate Specialization Process of Dermal Fibroblasts

Through UMAP dimensionality reduction cluster analysis, the gene expression profiles of different types of cells in yak hair follicle cycle changes were described in detail, and a series of cell type-specific genes were obtained. The results show that four different clusters were obtained via UMAP clustering for fibroblasts, which were clusters 6, 13, 14, and 16. Further GO enrichment analysis showed that the enriched GO pathways of clusters 13 and 14 were very similar, including response to stimulus, developmental process, and locomotion, but immune system process and rhythmic process were not enriched. For cluster 6, the characteristic genes expressed in these cells were significantly enriched in positive regulation of biological process, response to stimulus, regulation of biological process, metabolic process, cellular process, localization, and negative regulation of biological process, which shows the response of fibroblasts to external stimuli and the regulation of cell cycle, indicating a high degree of differentiation and tissue specificity (Figure 3a). At the same time, the results of the Circos map also show that there are many of the same genes and coenriched pathways between these different clusters, which also suggests that they are similar (Figure 3b). Based on the continuity of cells in the clustering results, we constructed pseudotime trajectories for four fibroblast-related cell populations to analyze the dynamic expression of genes during the differentiation of specific cell types. Using quasitime ordination analysis, one branch and three different states were obtained according to the gene expression pattern (Figure 3c), in which PF and RF were mainly concentrated in state 2, and myofibroblasts were mainly concentrated in state 3 (Figure 3d). The results show that the density of PF was high in the early stage and then gradually decreased, while the density of DF and myofibroblast was higher in the later stage (Figure 3e).
In addition, based on the quasitime sorting analysis, we analyzed the expression of four marker genes in the differentiation trajectory of fibroblasts, and found that the expression levels of TPM1 and NFIC were higher than those of FBLN2 during the differentiation process (Figure 4a). To further compare the changes in gene expression in different cell branches in different states, a series of genes were selected for differential analysis. For genes such as LRIG1 and TPM1, which are highly expressed before branching, the expression level is high in state 2 and low in state 1. The expression level of NFIC was upregulated in the process of cell specialization. FBLN2 was mainly concentrated in state 3, and its expression level was significantly lower than that of LRIG1, TPM1, and NFIC. We constructed a dynamic differentiation heatmap of differential genes (Figure 4b–d) and obtained four different gene sets according to their expression patterns. At the beginning of differentiation, it was mainly the specialization process of PF and RF cells. According to the pseudotemporal differentiation heatmap, CA6, KRT23, SLC15A1, KLK6, LY6E, CCHCR1, and GSDMA were highly expressed at this stage. For VIM, RGS1, CD74, and S100A4, the expression level was first upregulated and then downregulated. According to the gene expression matrix in the cluster heatmap results, the genes were divided into different genomes according to the expression pattern. Gene sets 1, 2, and 4 were mainly enriched in the cellular process, regulation of biological process, metabolic process, negative regulation of biological process, developmental regulation, locomotion, biological process involved in interspecies interaction between organisms, positive regulation of biological process, and response to stimulus (Figure 4e). In this study, the skin tissue of yaks was stained and analyzed via immunofluorescence (Figure 4f). The results show that VIM-positive cells were mainly concentrated in fibroblasts around hair follicles.

3.4. The Fate Specialization Process of Dermal Papilla Cells

DP cells were extracted from Seurat, and the Monocle algorithm was used to analyze the pseudotime differentiation trajectory of DP cells. We divided the specialization process into three branches. We analyzed the expression changes of DP marker genes in each state (Figure 5a). The expression levels of LEF1, MSX1, and HOXC13 in DP cells were first downregulated, then upregulated, and were highly expressed in state 6, among which HOXC13 was low in state 2, 3, and 4 (Figure 5b,c). Subsequently, we divided the branch-1-specific differential genes into four gene sets and performed expression analysis to further analyze the mechanism of DP cells in the hair follicle cycle transition process and the functional enrichment in the intermediate cell specialization process. The results show that gene set 1 was enriched in the negative regulation of external stimulus response, glycerolipid metabolism, epithelial cell proliferation, and cell–cell junction. Gene set 2 was enriched in intermediate filaments, intermediate filament cytoskeleton, low density lipoprotein particles, and chylomicrons. Gene set 3 was enriched in positive regulation of nociceptive response, cell–substrate junction assembly, skin development and integrin-mediated signaling pathway. Gene set 4 was enriched in hair follicle morphogenesis, positive regulation of interferon-β production, epidermal morphogenesis, and mitotic spindle assembly (Figure 5d).
We analyzed the branch-2-specific differential genes (Figure 6a). The results of enrichment analysis show that the second gene set was mainly enriched in lipid transfer activity, protein phosphatase inhibitor activity, and extracellular matrix binding. GO enrichment of the third gene set shows that it was mainly concentrated in pathways related to hair follicle development, such as epidermal development, extracellular matrix tissue, skin development, keratinocyte differentiation, hair cycle, epidermal cell differentiation, positive regulation of epithelial cell proliferation, and regulation of angiogenesis (Figure 6b). Among them, genes such as KRT17, SPINK5, and SOX9 are related to the cycle of the hair follicle cycle, and genes such as KRT17 and SOX9 are related to cell proliferation and hair follicle cycle development in the growth period. The high expression of these genes is more conducive to the regeneration of anagen hair follicles and the maintenance of cell homeostasis in degenerative hair follicles [47]. The fourth gene set was mainly enriched in the extrachromosomal kinetochores, DNA binding bending, and insulin-like growth factor binding, while IGFBP7 was involved in the related pathways of insulin-like growth factor binding. As a polypeptide growth factor, insulin-like growth factor plays an important role in cell anti-apoptosis, cell mitosis, and hair follicle dermal cell proliferation and differentiation [48]. Lee et al. [49] found that when follicular-keratinocyte-conditioned media (FKCM) was supplemented as a medium, proteins such as IGFBP7 had the function of enhancing hair induction, which enhanced the hair production of rat vibrissa DP cells and activated β-catenin and BMP signaling pathways. IGFBP3 is expressed in the dermal papilla, and its mRNA expression level is significantly increased in the early stages of degeneration and telogen compared to anagen [50]. We obtained data on different periods of the yak hair follicle cycle development from NCBI: PRJNA550233. The expression of IGFBP3 and IGFBP7 in different periods of the yak hair follicle cycle was observed via joint analysis. The results show that the fragments per kilobase of exon model per million fragments mapped (FPKM) values of IGFBP3 and IGFBP7 in the catagen and telogen groups were higher than those in the anagen group (Figure 6c). By further analyzing the expression of IGFBP3 and IGFBP7 in different states, it was found that the expression levels of IGFBP3 and IGFBP7 in dermal papilla cells were upregulated first and then downregulated, and the expression level was the lowest in state 6 (Figure 6d), which further indicates that the expression level of IGFBP3 and other genes in telogen were higher than those in anagen, which is consistent with the reported conclusions.
We analyzed the expression of branch-3-specific differential genes (Figure 7). The results show that the first gene set was mainly enriched in intermediate filament cytoskeleton, laminin binding, protein–lipid complex binding and extracellular matrix binding. The second gene set was mainly enriched in protein–DNA complex assembly, cell component decomposition involved in the execution phase of apoptosis, establishment of chromosome localization, and negative regulation of mid/late transition of cell cycle. Some researchers have analyzed the expression patterns of hair follicle cycle genes to show the differences in DP cells, and 825 differentially expressed genes were obtained from DP cells in the anagen phase and the telogen phase. Differential genes have functions such as cell cycle control, chromosome distribution, and cell division. It was found that secondary hair follicles in the middle stage of hair growth can not only induce the deformation of villi and the stagnation of villi growth, but also resist the apoptosis of hair follicles themselves and prepare for hair growth [51,52]. The third gene set was mainly enriched in the negative regulation of cell migration, the inflammatory response to wounds, the positive regulation of transforming growth factor β receptor signaling pathway, and the formation of animal organs. The fourth gene set was mainly enriched in items such as skin development, positive regulation of epithelial cell proliferation, keratinocyte differentiation, and hair cycle. The GO enrichment results of the fourth gene set were similar to those of the third gene set of the second branch.

3.5. CellChat Analysis

To understand cell interactions in the dermal fibroblast lineage, Cellchat [52] was used to infer the number and strength of interactions between different cell populations during the telogen and anagen of yak hair follicles. The results show that most of the interactions occurred between fibroblasts, IFE-DC, and IRS (Figure 8a). Then, the cell interaction between different cell types of fibroblasts was explored by the expression of ligand–receptor pairs. A cell communication analysis showed that RF has a high number and intensity of interactions with other cell types. Among them, the number of cell interactions between RF and PF and DF was greater, but the interaction intensity between RF and DF, DP, and myofibroblasts was significantly higher than that between RF and PF (Figure 8b). CellChat found multiple signaling pathways in the ligand–receptor pair of the fibroblast lineage, including the ncWNT, SPP1, and GRN pathways. In addition, the ncWNT signaling pathway from RF cells to DF cells was enriched, and the GRN signaling pathway was most abundant from myofibroblasts to PF cells (Figure 8c). Further analysis of the interaction of ligand–receptor pairs found that other cell types act on RF cells mainly through GRN-SORT1 ligand–receptor pairs. RF cells interact with other cell types mainly through SPP1-CD44 and GRN-SORT1 for cell communication (Figure 8d). Previous CellChat analysis of the ncWNT signaling network showed that ligands (WNT5A) and a fibroblast population mainly drive fibroblast-to-fibroblast, fibroblast-to-endothelial, and a small amount of fibroblast-to-medullary signals [53]. To reveal the effect of signaling pathways on intercellular communication, the interaction intensity of the three signaling pathways was demonstrated. Compared with other cells, RF cells were the most important cells for SPP1 signaling pathway output. In addition, the study also observed that the expression of the SPP1 and WNT11 genes was relatively high in the RF subgroup, and noncanonical WNT signals such as WNT11 were highly expressed in RF, suggesting that RF cells may communicate with DF cells through noncanonical WNT signaling pathways (Figure 8e,f).
In this study, immunofluorescence was used to analyze yak skin (Figure 9). The results show that MSX1 positive cells were expressed in the tissues of the hair bulb, especially in the lower DP cells. The expression level of DLX3 positive cells in DP cells was significantly lower than that of MSX1. Furthermore, we found that DLX3 also had a certain level of expression in the root sheath.

4. Discussion

Hair follicles are micro-organs with high cell heterogeneity. Hair follicle development starts from the embryonic period and eventually develops into an epithelial cell lineage and a dermal cell lineage composed of a variety of cells, and the whole development process is regulated by a variety of regulatory factors. Understanding the regulatory mechanism of hair follicle development and periodic growth of hair follicles is of great significance for the breeding of villus animal varieties and the improvement of villus yield. Fibroblasts are the most abundant cell type in the dermis of the skin. Thompson et al. [19] obtained 10 and 8 fibroblast clusters from scRNA-seq and scATAC-seq datasets, respectively, revealing that reticular fibroblasts are expected to follow the adipogenesis process, while neonatal fibroblasts may differentiate into the dermis. At the same time, during embryonic development, fibroblasts proliferate to fill the dermal structure and form dermal condensates along the regeneration trajectory [54,55,56]. Salzer et al. [57] found that due to development and aging, the skin maturation process may reduce the plasticity of the fibroblast population, and they found that upper dermal fibroblasts may easily establish cell contact and conduct signal transduction with adjacent cells, while lower dermal fibroblasts express more transcripts related to the extracellular matrix. Kim et al. [58] revealed the postnatal maturation trajectory of PFs via single-cell transcriptomics, and obtained cell clusters of the fibroblast lineage, including PF, RF, DP, and the dermal sheath, and two different differentiation trajectories in the upper fibroblast lineage were defined: Fb1 → Fb2 → Fb3 (FB maturation pathway) and Fb1 → Fb2 → Dp1 → Dp2 (DP regeneration pathway). There is a certain correlation between fibroblasts and hair follicles located in the epithelium and dermis in skin with rich hair [59]. PF has the special signal characteristics needed for hair follicle morphogenesis and coordination of hair growth [60]. At present, the internal mechanism of the spontaneous loss of regenerative capacity of dermal fibroblasts is still largely unknown. The differentiation of fibroblasts during the hair follicle cycle still needs further study.
At the beginning of telogen, there was no obvious expansion, apoptosis, or differentiation of hair follicles, and DP cells migrated to the bulge area near the hair follicle stem cells [61,62]. In telogen, the hair shaft is mature, and the distance between the DP and the hair bulb is shortened to the maximum extent. The transition from telogen to anagen is considered to be the activation of stem cells in hair follicle enlargement, the activation of dermal papilla cells, the interaction between the epidermis and the dermis, and the construction of the hair bulb, and the regeneration of secondary hair buds is the beginning of a new round of hair follicle cycle [63,64]. DP cells stop growing in the bulge area, so that DP cells and stem cells can interact directly. DP cells are essential for the regulation of stem cells and the initiation of the hair cycle. Studies have found that when the stem cell activator reaches a specific threshold, the growth period is initiated [5].
DP cells are the key cell type of hair follicle cycle transition, as they can coordinate mesenchymal–epithelial interactions and regulate hair follicle development and cycle transition [65]. Therefore, HFSCs and DP cells are necessary to initiate the first step of the new hair cycle [5,66]. The DP acts as the physical niche of progenitor cells during anagen and telogen, provides an opportunity to guide the proliferation of progenitor cells and the differentiation of their derivatives by inducing signal transduction [67]. Sun et al. [68] proposed that the DP sends signals to stem cells in the bulge region to initiate a new anagen. The DP expresses Wnts, R-spondins, FGFs, and Noggin, which can promote hair follicle growth and help initiate hair follicle regeneration [69,70]. DP cells can induce the transition from hair follicle anagen to catagen phase, and hair follicles in DP cells lacking β-catenin will enter the catagen in advance [71]. At present, the cycle transition at different stages of hair follicle development requires fine regulation between cell types such as DP. Therefore, if the molecular characteristics of DP cells during the transition of yak hair follicle cycle can be revealed in detail, it is of great significance to better understand the regulatory mechanism of the hair follicle cycle. In this study, scRNA-seq technology was used to analyze the molecular characteristics of DP cells in yak hair follicles from telogen to anagen. The expression profiles of three branch points were obtained, and the functional description of each lineage was described in detail. These findings provide new insights into hair regeneration.
During the transition of hair follicles from telogen to anagen, DP cells can secrete growth factors to participate in the regulation of hair follicle stem cell status and promote the cycle of hair follicles [72]. Chi et al. [73] proposed that the number of DP cells is directly related to the ability of hair follicles to initiate new hair growth: when the number of DP cells decreases below a specific threshold, hair follicles cannot initiate new hair cycles, but when a sufficient number of DPs is retained, hair follicles still have the ability to re-enter the growth phase. In mice, 20% of the hair follicles produce serrated hair in the primary integument, and more secondary hair types are produced in the second growth period, while more DP cells are needed for complex hair types, indicating that there may be functions related to hair type conversion in DP [21]. Hair follicle catagen is a dynamic transition period of the hair follicle cycle from anagen to telogen, which is regulated by molecules including FGF5 growth factor [74] and TGF-β pathway members [75,76]. IGF1 [77], PDGFB [78], WNT3A, WNT7A [79], and TAK1 [80] can maintain the anagen of hair follicles, of which only IGF1 [77] and PDGFA [81] are expressed in dermal papilla cells. There are two important signal transductions from the telogen to the anagen hair follicle cycle, in which β-catenin is the core pathway, and the switch that restarts the anagen hair follicle is the transient activation of β-catenin signaling [82,83,84,85]. It has also been found that RSPO1, similar to β-catenin, can also activate hair follicles into the anagen, but its mechanism is related to the activation of hair follicle stem cells [43]. BMP signaling plays an important role in the development of DP and the regulation of signal transmission [86]. In addition, it plays a role in inhibiting stem cells in the carina area [63]. BMP signaling can inhibit the transformation of hair follicles from telogen to anagen [87]. We found that the content of BMP2 was higher in the telogen of yak hair follicles. It is worth discussing that whether there are similarities and differences in the mechanism of action of BMP and other signals in yak hair follicles at different ages and between yak and other model animals is also a direction for further research in the future. We analyzed the differential gene expression of DP cells, and found that these genes were differentially expressed in different DP cell cycles, and reflected the differences in time and space, which provided an important reference for future research.
From the GO enrichment entries, it can be seen that DP cells are enriched in pathways such as the cell–cell junction, laminin binding, protein–DNA complex assembly, and protein–lipid complex, indicating that dermal papilla cells are responsible for the expression and transmission of genetic information. By analyzing the different branches, it was found that the first branch of dermal papilla cells mainly focused on cell proliferation and the function of participating in cell junction information transmission, which also reflected the importance of DP cells as the control center of hair follicle growth. DP cells can subtly control the continuous growth of mammalian hair by releasing a large number of molecular signals [88]. During the whole hair cycle, the DP is in a favorable position, where it can provide induction signals and guide the cell activity of hair follicles. The second branch focuses on the pathways of cell differentiation, hair cycle and angiogenesis. These pathways may be involved in the downward migration of hair follicles from telogen to anagen to the dermis, which is connected with capillaries and provides nutrients for hair growth. In the third branch, not only were cell connection information transmission and other related functions enriched, but apoptosis-related functions were also enriched in the second gene set. Yang et al. [32] found that the second branch of cashmere goat DP cells was enriched in the apoptosis process via a pseudotemporal analysis of cashmere goat DP cells. They speculated that the regulation of normal gene transcription into pseudogenes in the process of apoptosis coordinated hair follicle apoptosis. We analyzed the differential gene expression and potential functions of DP cells, but more details, such as the spatial distribution pattern of gene expression in DP cells, have not been studied and described in detail in mice and yaks. Future research may focus on such a theme, which may provide new insights for our subsequent research on hair follicle biology and hair follicle regeneration.
A CellChat analysis revealed a high number and intensity of interactions between RF and other cell types in the dermal fibroblast lineage. Notably, RF was identified as the most important cell for SPP1 signaling pathway output. This means that SPP1 plays an important role in fibroblasts and shows high expression characteristics under matrix activation, which further promotes the high expression of fibrosis genes [89]. By using Tyr-NrasQ61K; Spp1−/−mice, researchers found that SPP1 loss-of-function mutations are sufficient to reverse the static hair cycle in the skin of congenital nevus and that SPP1 can induce new hair growth [90]. Periodic hair growth and dormancy are complex processes regulated by a variety of internal and external factors, and SPP1 is a key mediator, which indicates that high expression of SPP1 in fibroblasts may play an important role in regulating hair cycle and wound repair. Liaw et al. [91] confirmed the role of SPP1 in the process of wound repair; by simultaneously knocking out the SPP1 gene in the mouse skin incision model, it was observed that the matrix disintegration and collagen fiber diameter decreased and the incision could not be repaired, suggesting the importance of OPN in the process of wound repair. In addition, inflammatory-cell-mediated signaling can trigger the upregulation of SPP1 in wound fibroblasts, which may delay the repair process and may also partially lead to fibrosis caused by wound healing [92].

5. Conclusions

In summary, this study used scRNA-seq technology to draw a cell transcription map in the skin, systematically described the changes in the expression of time-like genes during dermal fibroblast lineage and DP cell specialization, constructed an intercellular communication network, and demonstrated the interaction between identified cell types. The current research also highlights the previously neglected cell fate determination process in the cycle of yak hair follicles. In addition, by dissecting the heterogeneity of fibroblast subtypes, our study highlights that different fibroblast signals coordinate the asynchronous development of different cell types. The results of this study can be used to promote further research on the molecular regulatory mechanism of hair follicles and provide new insights into the functional analysis of dermal cell lineage during the transition of yak hair follicles from telogen to anagen.

Author Contributions

Conceptualization, P.Y., H.P. and Q.Z.; methodology, Q.Z.; software, Q.Z. and N.Y.; validation, P.Y., H.P. and Q.Z.; formal analysis, Q.Z., T.W. and N.Y.; investigation, Q.Z., X.W., M.C. and C.M.; resources, N.Y., P.B., S.K. and C.L.; data curation, M.C., X.G. and P.B.; writing—original draft preparation, Q.Z., T.W. and C.M.; writing—review and editing, X.W., X.G., S.K. and C.L.; visualization, H.P. and P.Y.; supervision, P.Y.; project administration, P.Y.; funding acquisition, P.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research Program (2021YFD1600201), the National Natural Science Foundation of China (31972561), the National Beef Cattle Industry Technology and System (CARS-37), the Agricultural Science and Technology Innovation Program (25-LZIHPS-01), the Science and Technology Aid Qinghai Cooperation Special Project (2020-QY-212), the Shenzhen Science and Technology Program (Grant No. KCXFZ20201221173205015), and the Science and Technology program of Gansu Province (20JR5RA580).

Institutional Review Board Statement

All skin samples were collected from anesthetized yaks. All procedures involving animals were performed according to the guidelines of the China Council on Animal Care and the Ministry of Agriculture of the People’s Republic of China, as well as the Animal Care and Use Committee of the Lanzhou Institute of Husbandry and Pharmaceutical Sciences Chinese Academy of Agricultural Sciences (No. LIHPS-CAAS- 2017-115, date of approval 2017).

Informed Consent Statement

Not applicable.

Data Availability Statement

All the sequencing data used in this research are deposited in NCBI GEO databases under accession number: GSE205456. Data on telogen are scheduled to be released on 30 May 2024.

Acknowledgments

We thank the Berry Genomics Corporation (Beijing, China) for its help in sequencing the samples. We thank our laboratory colleagues, Yi Wu and Qi Bao, for assistance with the data analysis.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bao, P.; Luo, J.; Liu, Y.; Chu, M.; Ren, Q.; Guo, X.; Tang, B.; Ding, X.; Qiu, Q.; Pan, H.; et al. The seasonal development dynamics of the yak hair cycle transcriptome. BMC Genom. 2020, 21, 355. [Google Scholar] [CrossRef]
  2. Milner, Y.; Sudnik, J.; Filippi, M.; Kizoulis, M.; Kashgarian, M.; Stenn, K. Exogen, shedding phase of the hair growth cycle: Characterization of a mouse model. J. Investig. Dermatol. 2002, 119, 639–644. [Google Scholar] [CrossRef]
  3. Alonso, L.; Fuchs, E. The hair cycle. J. Cell Sci. 2006, 119, 391–393. [Google Scholar] [CrossRef]
  4. Paus, R.; Foitzik, K. In search of the “hair cycle clock”: A guided tour. Differentiation 2004, 72, 489–511. [Google Scholar] [CrossRef]
  5. Schneider, M.; Schmidt-Ullrich, R.; Paus, R. The hair follicle as a dynamic miniorgan. Curr. Biol. 2009, 19, 132–142. [Google Scholar] [CrossRef]
  6. Lee, J.; Tudorita, T. Hairy tale of signaling in hair follicle development and cycling. Semin. Cell Dev. Biol. 2012, 23, 906–916. [Google Scholar]
  7. Paus, R.; Cotsarelis, G. The biology of hair follicles. N. Engl. J. Med. 1999, 341, 491–497. [Google Scholar] [CrossRef]
  8. Botchkarev, V.; Kishimoto, J. Molecular control of epithelial-mesenchymal interactions during hair follicle cycling. J. Investig. Dermatol. Symp. Proc. 2003, 8, 46–55. [Google Scholar] [CrossRef]
  9. Jorgenson, A.J.; Choi, K.M.; Sicard, D. TAZ activation drives fibroblast spheroid growth, expression of profibrotic paracrine signals, and context-dependent ECM gene expression. Am. J. Physiol. Cell Physiol. 2017, 312, 277–285. [Google Scholar] [CrossRef]
  10. Sorrell, J.M.; Caplan, A.I. Fibroblast heterogeneity: More than skin deep. J. Cell Sci. 2004, 117, 667–675. [Google Scholar] [CrossRef]
  11. Driskell, R.R.; Lichtenberger, B.M.; Hoste, E.; Kretzschmar, K.; Simons, B.D.; Charalambous, M.; Ferron, S.R.; Herault, Y.; Pavlovic, G.; Ferguson-Smith, A.C.; et al. Distinct fibroblast lineages determine dermal architecture in skin development and repair. Nature 2013, 504, 277–281. [Google Scholar] [CrossRef]
  12. Mine, S.; Fortunel, N.O.; Pageon, H.; Asselineau, D. Aging alters functionally human dermal papillary fibroblasts but not reticular fibroblasts: A new view of skin morphogenesis and aging. PLoS ONE 2008, 3, e4066. [Google Scholar] [CrossRef]
  13. Janson, D.G.; Saintigny, G.; Adrichem, A.V.; Mahé, C.; Ghalbzouri, A.E. Different gene expression patterns in human papillary and reticular fibroblasts. J. Investig. Dermatol. 2012, 132, 2565–2572. [Google Scholar] [CrossRef]
  14. Lewis, D.A.; Travers, J.B.; Machado, C.; Somani, A.; Spandau, D.F. Reversing the aging stromal phenotype prevents carcinoma initiation. Aging 2011, 3, 407–416. [Google Scholar] [CrossRef]
  15. Mascharak, S.; Desjardins-Park, H.E.; Davitt, M.F.; Griffin, M.; Borrelli, M.R.; Moore, A.L.; Chen, K.; Duoto, B.; Chinta, M.; Foster, D.S.; et al. Preventing engrailed-1 activation in fibroblasts yields wound regeneration without scarring. Science 2021, 372, eaba2374. [Google Scholar] [CrossRef]
  16. Plikus, M.V.; Guerrero-Juarez, C.F.; Ito, M. Regeneration of fat cells from myofibroblasts during wound healing. Science 2017, 355, 748–752. [Google Scholar] [CrossRef]
  17. Guerrero-Juarez, C.F.; Dedhia, P.H.; Jin, S.; Ruiz-Vega, R.; Ma, D.; Liu, Y.; Yamaga, K.; Shestova, O.; Gay, D.L.; Yang, Z.; et al. Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds. Nat. Commun. 2019, 10, 18–32. [Google Scholar] [CrossRef]
  18. Vorstandlechner, V.; Laggner, M.; Kalinina, P.; Haslik, W.; Radtke, C.; Shaw, L.; Lichtenberger, B.M.; Tschachler, E.; Ankersmit, H.J.; Mildner, M. Deciphering the functional heterogeneity of skin fibroblasts using single-cell RNA sequencing. FASEB J. 2020, 34, 3677–3692. [Google Scholar] [CrossRef]
  19. Thompson, S.M.; Phan, Q.M.; Winuthayanon, S.; Driskell, I.M.; Driskell, R.R. Parallel single-cell multiomics analysis of neonatal skin reveals the transitional fibroblast states that restrict differentiation into distinct fates. J. Investig. Dermatol. 2022, 142, 1812–1823. [Google Scholar] [CrossRef]
  20. Lu, K.; Han, Q.; Ma, Z.; Yan, Q.; Pei, Y.; Shi, P.; Zhang, J.; Rong, K.; Ma, K.; Li, P.; et al. Injectable platelet rich fibrin facilitates hair follicle regeneration by promoting human dermal papilla cell proliferation, migration, and trichogenic inductivity. Exp. Cell Res. 2021, 409, 112888. [Google Scholar] [CrossRef]
  21. Chi, W.; Wu, E.; Morgan, B.A. Dermal papilla cell number specifies hair size, shape and cycling and its reduction causes follicular decline. Development 2013, 140, 1676–1683. [Google Scholar] [CrossRef]
  22. Tan, J.J.Y.; Common, J.E.; Wu, C. Keratinocytes maintain compartmentalization between dermal papilla and fibroblasts in 3D heterotypic tri-cultures. Cell Prolif. 2019, 52, e12668. [Google Scholar] [CrossRef] [PubMed]
  23. Rompolas, P.; Deschene, E.R.; Zito, G.; Gonzalez, D.G.; Saotome, I.; Haberman, A.M.; Greco, V. Live imaging of stem cell and progeny behaviour in physiological hair-follicle regeneration. Nature 2012, 487, 496–499. [Google Scholar] [CrossRef] [PubMed]
  24. La, Y.; Ma, X.; Bao, P.; Chu, M.; Yan, P.; Guo, X.; Liang, C. Identification and characterization of piwi-interacting RNAs for early testicular development in yak. Int. J. Mol. Sci. 2022, 23, 12320. [Google Scholar] [CrossRef] [PubMed]
  25. Zhang, X.; Bao, P.; Ye, N.; Zhou, X.; Zhang, Y.; Liang, C.; Guo, X.; Chu, M.; Pei, J.; Yan, P. Identification of the key genes associated with the yak hair follicle cycle. Genes 2021, 13, 32. [Google Scholar] [CrossRef] [PubMed]
  26. Zheng, Q.; Ye, N.; Bao, P.; Zhang, X.; Wang, F.; Ma, L.; Chu, M.; Guo, X.; Liang, C.; Pan, H.; et al. Construction of transcriptome atlas of white yak hair follicle during anagen and catagen using single-cell RNA sequencing. BMC Genom. 2022, 23, 813. [Google Scholar] [CrossRef] [PubMed]
  27. He, W.; Ye, J.; Xu, H.; Lin, Y.; Zheng, Y. Differential expression of α6 and β1 integrins reveals epidermal heterogeneity at single-cell resolution. J. Cell Biochem. 2020, 121, 2664–2676. [Google Scholar] [CrossRef]
  28. Joost, S.; Zeisel, A.; Jacob, T.; Sun, X.; Manno, G.L.; Lönnerberg, P.; Linnarsson, S.; Kasper, M. Single-cell transcriptomics reveals that differentiation and spatial signatures shape epidermal and hair follicle heterogeneity. Cell Syst. 2016, 3, 221–237. [Google Scholar] [CrossRef]
  29. Schlake, T. Segmental Igfbp5 expression is specifically associated with the bent structure of zigzag hairs. Mech. Dev. 2005, 122, 988–997. [Google Scholar] [CrossRef]
  30. Yang, H.; Adam, R.C.; Ge, Y.; Hua, Z.L.; Fuchs, E. Epithelial-mesenchymal micro-niches govern stem cell lineage choices. Cell 2017, 169, 483–496. [Google Scholar] [CrossRef]
  31. Greco, V.; Chen, T.; Rendl, M.; Schober, M.; Pasolli, H.A.; Stokes, N.; Cruz-Racelis, J.D.; Fuchs, E. A two-step mechanism for stem cell activation during hair regeneration. Cell Stem Cell 2009, 4, 155–169. [Google Scholar] [CrossRef] [PubMed]
  32. Yang, F.; Li, R.; Zhao, C.; Che, T.; Guo, J.; Xie, Y.; Wang, Z.; Li, J.; Liu, Z. Single-cell sequencing reveals the new existence form of dermal papilla cells in the hair follicle regeneration of cashmere goats. Genomics 2022, 114, 110316. [Google Scholar] [CrossRef] [PubMed]
  33. Gu, L.H.; Coulombe, P.A. Keratin function in skin epithelia: A broadening palette with surprising shades. Curr. Opin. Cell Biol. 2007, 19, 13–23. [Google Scholar] [CrossRef] [PubMed]
  34. Ge, W.; Zhang, W.; Zhang, Y.; Zheng, Y.; Li, F.; Wang, S.; Liu, J.; Tan, S.; Yan, Z.; Wang, L.; et al. A single-cell transcriptome atlas of cashmere goat hair follicle morphogenesis. Genom. Proteom. Bioinform. 2021, 19, 437–451. [Google Scholar] [CrossRef] [PubMed]
  35. Rezza, A.; Wang, Z.; Sennett, R.; Qiao, W.; Wang, D.; Heitman, N.; Mok, K.W.; Clavel, C.; Yi, R.; Zandstra, P. Signaling networks among stem cell precursors, transit-amplifying progenitors, and their niche in developing hair follicles. Cell Rep. 2016, 14, 3001–3018. [Google Scholar] [CrossRef] [PubMed]
  36. Morita, R.; Sanzen, N.; Sasaki, H.; Hayashi, T.; Umeda, M.; Yoshimura, M.; Yamamoto, T. Tracing the origin of hair follicle stem cells. Nature 2021, 594, 547–552. [Google Scholar] [CrossRef] [PubMed]
  37. Chovatiya, G.; Ghuwalewala, S.; Walter, L.D.; Cosgrove, B.D.; Tumbar, T. High resolution single cell transcriptomics reveals heterogeneity of self-renewing hair follicle stem cells. Exp. Dermatol. 2020, 30, 457–471. [Google Scholar] [CrossRef]
  38. Wu, S.; Yu, Y.; Liu, C.; Zhang, X.; Zhu, P.; Peng, Y.; Yan, X. Single-cell transcriptomics reveals lineage trajectory of human scalp hair follicle and informs mechanisms of hair graying. Cell Discov. 2022, 8, 49. [Google Scholar] [CrossRef]
  39. Joost, S.; Annusver, K.; Jacob, T.; Sun, X.; Dalessandri, T.; Sivan, U.; Sequeira, I. The molecular anatomy of mouse skin during hair growth and rest. Cell Stem Cell 2020, 26, 441–457. [Google Scholar] [CrossRef]
  40. Morgan, H.J.; Benketah, A.; Olivero, C.; Rees, E.; Ziaj, S.; Mukhtar, A.; Lanfredini, S.; Patel, G.K. Hair follicle differentiation-specific keratin expression in human basal cell carcinoma. Clin. Exp. Dermatol. 2020, 45, 417–425. [Google Scholar] [CrossRef]
  41. Sano, S.; Kira, M.; Takagi, S.; Yoshikawa, K.; Takeda, J.; Itami, S. Two distinct signaling pathways in hair cycle induction: Stat3-dependent and -independent pathways. Proc. Natl. Acad. Sci. USA 2000, 97, 13824–13829. [Google Scholar] [CrossRef]
  42. Huelsken, J.; Vogel, R.; Erdmann, B.; Cotsarelis, G.; Birchmeier, W. beta-Catenin controls hair follicle morphogenesis and stem cell differentiation in the skin. Cell 2001, 105, 533–545. [Google Scholar] [CrossRef] [PubMed]
  43. Li, N.; Liu, S.; Zhang, H.; Deng, Z.; Zhao, H.; Zhao, Q.; Lei, X. Exogenous R-spondin1 induces precocious telogen-to-anagen transition in mouse hair follicles. Int. J. Mol. Sci. 2016, 17, 582. [Google Scholar] [CrossRef]
  44. Chua, A.W.C.; Ma, D.; Gan, S.U.; Fu, Z.; Han, H.C.; Song, C.; Sabapathy, K.; Phan, T.T. The role of R-spondin2 in keratinocyte proliferation and epidermal thickening in keloid scarring. J. Investig. Dermatol. 2011, 131, 644–654. [Google Scholar] [CrossRef]
  45. Parma, P.; Radi, O.; Vidal, V.; Chaboissier, M.C.; Dellambra, E.; Valentini, S.; Guerra, L.; Schedl, A.; Camerino, G. R-spondin1 is essential in sex determination, skin differentiation and malignancy. Nat. Genet. 2006, 38, 1304–1309. [Google Scholar] [CrossRef]
  46. Guo, H.; Yang, K.; Deng, F.; Ye, J.; Xing, Y.; Li, Y.; Lian, X.; Yang, T. Wnt3a promotes melanin synthesis of mouse hair follicle melanocytes. Biochem. Biophys. Res. Commun. 2012, 420, 799–804. [Google Scholar] [CrossRef] [PubMed]
  47. Wang, J. The Role of H3K27me3 on the Cyclic Regenesis of Hair Follicle from Shaanbei White Cashmere Goats; Northwest A& F University: Xianyang, China, 2021. [Google Scholar]
  48. Wudubala, L.Y.; Wu, T.; Ma, Y.; Gaowa; Zhao, X.; Yu, P.; Liu, B. Research progress on regulative role of IGF system in cashmere growth and reproductive activity of cashmere goat. Anim. Husb. Feed Sci. 2020, 41, 28–33. [Google Scholar]
  49. Lee, M.H.; Im, S.; Shin, S.H.; Kwack, M.H.; Jun, S.; Kim, M.K.; Kim, J.C.; Sung, Y.K. Conditioned media obtained from human outer root sheath follicular keratinocyte culture activates signalling pathways that contribute to maintenance of hair-inducing capacity and increases trichogenicity of cultured dermal cells. Exp. Dermatol. 2012, 21, 793–795. [Google Scholar] [CrossRef]
  50. Schlake, T.; Beibel, M.; Weger, N.; Boehm, T. Major shifts in genomic activity accompany progression through different stages of the hair cycle. Gene Expr. Patterns 2004, 4, 141–152. [Google Scholar] [CrossRef]
  51. Zhu, B.; Xu, T.; Yuan, J.; Guo, X.; Liu, D. Transcriptome sequencing reveals differences between primary and secondary hair follicle-derived dermal papilla cells of the cashmere goat (capra hircus). Physiol. Genom. 2014, 46, 104–111. [Google Scholar] [CrossRef]
  52. Yang, F. Study on the Action Mechanism of Hair Follicles Cells in Inner Mongolia Cashmere Goat; Inner Mongolia Agricultural University: Hohhot, China, 2020. [Google Scholar]
  53. Jin, S.; Guerrero, J.C.F.; Zhang, L. Inference and analysis of cell-cell communication using CellChat. Nat. Commun. 2021, 12, 1088. [Google Scholar] [CrossRef] [PubMed]
  54. Chen, D.; Jarrell, A.; Guo, C.; Lang, R.; Atit, R. Dermal β-catenin activity in response to epidermal Wnt ligands is required for fibroblast proliferation and hair follicle initiation. Development 2012, 139, 1522–1533. [Google Scholar] [CrossRef] [PubMed]
  55. Gupta, K.; Levinsohn, J.; Linderman, G.; Chen, D.; Sun, T.Y.; Dong, D.; Taketo, M.M. Single-cell analysis reveals a hair follicle dermal niche molecular differentiation trajectory that begins prior to morphogenesis. Dev. Cell 2018, 48, 17–31. [Google Scholar] [CrossRef]
  56. Mok, K.; Saxena, N.; Heitman, N.; Grisanti, L.; Srivastava, D.; Muraro, M.J.; Jacob, T. Dermal condensate niche fate specification occurs prior to formation and is placode progenitor dependent. Dev. Cell 2018, 48, 32–48. [Google Scholar] [CrossRef] [PubMed]
  57. Salzer, M.C.; Lafzi, A.; Berenguer-Llergo, A.; Youssif, C.; Castellanos, A.; Solanas, G.; Peixoto, F.O. Identity noise and adipogenic traits characterize dermal fibroblast aging. Cell 2018, 175, 1575–1590. [Google Scholar] [CrossRef]
  58. Kim, J.Y.; Park, M.; Ohn, J.; Seong, R.H.; Chung, J.H.; Kim, K.H.; Jo, S.J.; Kwon, O. Twist2-driven chromatin remodeling governs the postnatal maturation of dermal fibroblasts. Cell Rep. 2022, 39, 110821. [Google Scholar] [CrossRef] [PubMed]
  59. Ohyama, M.; Zheng, Y.; Paus, R.; Stenn, K.S. The mesenchymal component of hair follicle neogenesis: Background, methods and molecular characterization. Exp. Dermatol. 2010, 19, 89–99. [Google Scholar] [CrossRef]
  60. Sennett, R.; Rendl, M. Mesenchymal-epithelial interactions during hair follicle morphogenesis and cycling. Semin. Cell Dev. Biol. 2012, 23, 917–927. [Google Scholar] [CrossRef]
  61. Taylor, G.; Lehrer, M.S.; Jensen, P.J.; Sun, T.T.; Lavker, R.M. Involvement of follicular stem cells in forming not only the follicle but also the epidermis. Cell 2000, 102, 451–461. [Google Scholar] [CrossRef]
  62. Oshima, H.; Rochat, A.; Kedzia, C.; Kobayashi, K.; Barrandon, Y. Morphogenesis and renewal of hair follicles from adult multipotent stem cells. Cell 2001, 104, 233–245. [Google Scholar] [CrossRef]
  63. Blanpain, C.; Lowry, W.E.; Geoghegan, A.; Polak, L.; Fuchs, E. Self-renewal, multipotency, and the existence of two cell populations within an epithelial stem cell niche. Cell 2004, 118, 635–648. [Google Scholar] [CrossRef] [PubMed]
  64. Tumbar, T.; Guasch, G.; Greco, V.; Blanpain, C.; Lowry, W.E.; Rendl, M.; Fuchs, E. Defining the epithelial stem cell niche in skin. Science 2004, 303, 359–363. [Google Scholar] [CrossRef]
  65. Namekata, M.; Yamamoto, M.; Goitsuka, R. Nuclear localization of meis1 in dermal papilla promotes hair matrix cell proliferation in the anagen phase of hair cycle. Biochem. Biophys. Res. Commun. 2019, 519, 727–733. [Google Scholar] [CrossRef] [PubMed]
  66. Yi, R.; Fuchs, E. MicroRNA-mediated control in the skin. Cell Death Differ. 2010, 17, 229–235. [Google Scholar] [CrossRef] [PubMed]
  67. Morgan, B.A. The dermal papilla: An instructive niche for epithelial stem and progenitor cells in development and regeneration of the hair follicle. Cold Spring Harb. Perspect. Med. 2014, 4, a015180. [Google Scholar] [CrossRef]
  68. Sun, T.T.; Cotsarelis, G.; Lavker, R.M. Hair follicular stem cells: The bulge-activation hypothesis. J. Investig. Dermatol. 1991, 96, 77–78. [Google Scholar] [CrossRef] [PubMed]
  69. Reddy, S.; Andl, T.; Bagasra, A.; Lu, M.M.; Epstein, D.J.; Morrisey, E.E.; Millar, S.E. Characterization of Wnt gene expression in developing and postnatal hair follicles and identification of Wnt5a as a target of Sonic hedgehog in hair follicle morphogenesis. Mech. Dev. 2001, 107, 69–82. [Google Scholar] [CrossRef]
  70. Rendl, M.; Lewis, L.; Fuchs, E. Molecular dissection of mesenchymal-epithelial interactions in the hair follicle. PLoS Biol. 2005, 3, e331. [Google Scholar] [CrossRef]
  71. Enshell-Seijffers, D.; Lindon, C.; Kashiwagi, M.; Morgan, B.A. beta-catenin activity in the dermal papilla regulates morphogenesis and regeneration of hair. Dev. Cell 2010, 18, 633–642. [Google Scholar] [CrossRef]
  72. Rahmani, W.; Abbasi, S.; Hagner, A.; Raharjo, E.; Kumar, R.; Hotta, A.; Magness, S. Hair follicle dermal stem cells regenerate the dermal sheath, repopulate the dermal papilla, and modulate hair type. Dev. Cell 2014, 31, 543–558. [Google Scholar] [CrossRef]
  73. Chi, W.Y.; Enshell-Seijffers, D.; Morgan, B.A. De novo production of dermal papilla cells during the anagen phase of the hair cycle. J. Investig. Dermatol. 2010, 130, 2664–2666. [Google Scholar] [CrossRef]
  74. Hébert, J.M.; Rosenquist, T.; Götz, J.; Martin, G.R. FGF5 as a regulator of the hair growth cycle: Evidence from targeted and spontaneous mutations. Cell 1994, 78, 1017–1025. [Google Scholar] [CrossRef]
  75. Foitzik, K.; Lindner, G.; Mueller-Roever, S.; Maurer, M.; Botchkareva, N.; Botchkarev, V.; Handjiski, B.; Metz, M. Control of murine hair follicle regression (catagen) by TGF-beta1 in vivo. FASEB J. 2000, 14, 752–760. [Google Scholar] [CrossRef] [PubMed]
  76. Andl, T.; Ahn, K.; Kairo, A.; Chu, E.Y.; Wine-Lee, L.; Reddy, S.T. Epithelial Bmpr1a regulates differentiation and proliferation in postnatal hair follicles and is essential for tooth development. Development 2004, 131, 2257–2268. [Google Scholar] [CrossRef] [PubMed]
  77. Ahn, S.; Pi, L.; Hwang, S.T.; Lee, W. Effect of IGF-I on hair growth is related to the anti-apoptotic effect of IGF-I and up-regulation of PDGF-A and PDGF-B. Ann. Dermatol. 2012, 24, 26–31. [Google Scholar] [CrossRef] [PubMed]
  78. Tomita, Y.; Akiyama, M.; Shimizu, H. PDGF isoforms induce and maintain anagen phase of murine hair follicles. J. Dermatol. Sci. 2006, 43, 105–115. [Google Scholar] [CrossRef] [PubMed]
  79. Kishimoto, J.; Burgeson, R.E.; Morgan, B.A. Wnt signaling maintains the hair-inducing activity of the dermal papilla. Genes Dev. 2000, 14, 1181–1185. [Google Scholar] [CrossRef]
  80. Sayama, K.; Kajiya, K.; Sugawara, K.; Sato, S.; Hirakawa, S.; Shirakata, Y.; Hanakawa, Y. Inflammatory mediator TAK1 regulates hair follicle morphogenesis and anagen induction shown by using keratinocyte-specific TAK1-deficient mice. PLoS ONE 2010, 5, e11275. [Google Scholar] [CrossRef]
  81. Kamp, H.; Geilen, C.C.; Sommer, C.; Blume-Peytavi, U. Regulation of PDGF and PDGF receptor in cultured dermal papilla cells and follicular keratinocytes of the human hair follicle. Exp. Dermatol. 2003, 12, 662–672. [Google Scholar] [CrossRef]
  82. Celso, C.L.; Prowse, D.M.; Watt, F.M. Transient activation of beta-catenin signalling in adult mouse epidermis is sufficient to induce new hair follicles but continuous activation is required to maintain hair follicle tumours. Development 2004, 131, 1787–1799. [Google Scholar] [CrossRef]
  83. Lowry, W.E.; Blanpain, C.; Nowak, J.A.; Guasch, G.; Lewis, L.; Fuchs, E. Defining the impact of beta-catenin/Tcf transactivation on epithelial stem cells. Genes Dev. 2005, 19, 1596–1611. [Google Scholar] [CrossRef]
  84. Zhang, Y.; Andl, T.; Yang, S.H.; Teta, M.; Liu, F.; Seykora, J.T.; Tobias, J.W. Activation of beta-catenin signaling programs embryonic epidermis to hair follicle fate. Development 2008, 135, 2161–2172. [Google Scholar] [CrossRef] [PubMed]
  85. Mater, D.V.; Kolligs, F.T.; Dlugosz, A.A.; Fearon, E.R. Transient activation of beta-catenin signaling in cutaneous keratinocytes is sufficient to trigger the active growth phase of the hair cycle in mice. Genes Dev. 2003, 17, 1219–1224. [Google Scholar] [CrossRef] [PubMed]
  86. Rendl, M.; Polak, L.; Fuchs, E. BMP signaling in dermal papilla cells is required for their hair follicle-inductive properties. Genes Dev. 2008, 22, 543–557. [Google Scholar] [CrossRef] [PubMed]
  87. Zhang, J.; He, X.C.; Tong, W.; Johnson, T.; Wiedemann, L.M.; Mishina, Y.; Feng, J.Q.; Li, L. Bone morphogenetic protein signaling inhibits hair follicle anagen induction by restricting epithelial stem/progenitor cell activation and expansion. Stem Cells 2006, 24, 2826–2839. [Google Scholar] [CrossRef] [PubMed]
  88. Ma, S.; Wang, Y.; Zhou, G.; Ding, Y.; Yang, Y.; Wang, X.; Zhang, E.; Chen, Y. Synchronous profiling and analysis of mRNAs and ncRNAs in the dermal papilla cells from cashmere goats. BMC Genom. 2019, 20, 512. [Google Scholar] [CrossRef]
  89. Cheng, J.; Wu, H.; Xie, C.; He, Y.; Mou, R.; Zhang, H.; Yang, Y.; Xu, Q. Single cell mapping of large and small arteries during hypertensive aging. J. Gerontol. A. Biol. Sci. Med. Sci. 2023, 188, glad188. [Google Scholar] [CrossRef]
  90. Wang, X.; Ramos, R.; Phan, A.Q.; Yamaga, K.; Flesher, J.L.; Jiang, S.; Oh, J.W.; Jin, S.; Jahid, S.; Kuan, C.H.; et al. Signalling by senescent melanocytes hyperactivates hair growth. Nature 2023, 618, 808–817. [Google Scholar] [CrossRef]
  91. Liaw, L.; Birk, D.E.; Ballas, C.B.; Whitsitt, J.S.; Davidson, J.M.; Hogan, B.L. Altered wound healing in mice lacking a functional osteopontin gene (spp1). J. Clin. Investig. 1998, 101, 1468–1478. [Google Scholar] [CrossRef]
  92. Mori, R.; Shaw, T.J.; Martin, P. Molecular mechanisms linking wound inflammation and fibrosis: Knockdown of osteopontin leads to rapid repair and reduced scarring. J. Exp. Med. 2008, 205, 43–51. [Google Scholar] [CrossRef]
Animals 13 03818 g001
Figure 1. Single-cell UMAP cluster analysis. (a) Data integration analysis and identification of major cell types; (b) main cell type proportions during telogen and anagen; (c) the expression of different types of cell-specific marker molecules in all cells; (d) the expression of characteristic genes in different cell types.
Figure 1. Single-cell UMAP cluster analysis. (a) Data integration analysis and identification of major cell types; (b) main cell type proportions during telogen and anagen; (c) the expression of different types of cell-specific marker molecules in all cells; (d) the expression of characteristic genes in different cell types.
Animals 13 03818 g001
Animals 13 03818 g002
Figure 2. The expression of key genes in anagen and telogen of hair follicles.
Figure 2. The expression of key genes in anagen and telogen of hair follicles.
Animals 13 03818 g002
Animals 13 03818 g003
Figure 3. Comparison and pseudotiming analysis of GO pathways enriched in different clusters of fibroblast lineage. (a) Comparison of GO pathways enriched in different clusters of fibroblast lineage; (b) The Circos map shows the GO terms shared by the four gene sets; (c) construction of the pseudotime differentiation trajectory of fibroblasts; (d) cell branch trajectory map; (e) fibroblast tree diagram and cell density diagram.
Figure 3. Comparison and pseudotiming analysis of GO pathways enriched in different clusters of fibroblast lineage. (a) Comparison of GO pathways enriched in different clusters of fibroblast lineage; (b) The Circos map shows the GO terms shared by the four gene sets; (c) construction of the pseudotime differentiation trajectory of fibroblasts; (d) cell branch trajectory map; (e) fibroblast tree diagram and cell density diagram.
Animals 13 03818 g003
Animals 13 03818 g004
Figure 4. Molecular characteristics in the process of fibroblast specialization. (a) Differential expression of characteristic genes in branches; (b) changes in the expression of characteristic genes over time; (c) the expression of characteristic genes in the lower state; (d) gene dynamic changes in the process of cell specialization; (e) comparison of GO terms in different gene sets; (f) immunofluorescence analysis.
Figure 4. Molecular characteristics in the process of fibroblast specialization. (a) Differential expression of characteristic genes in branches; (b) changes in the expression of characteristic genes over time; (c) the expression of characteristic genes in the lower state; (d) gene dynamic changes in the process of cell specialization; (e) comparison of GO terms in different gene sets; (f) immunofluorescence analysis.
Animals 13 03818 g004
Animals 13 03818 g005
Figure 5. DP cell differentiation trajectory and dynamic changes in the branch 1 gene. (a) Construction of the DP pseudotime differentiation trajectory; (b) changes in the expression of LEF1, MSX1, and HOXC13 in pseudotime; (c) heatmap displaying the branch 1 expression pattern; (d) GO enrichment analysis of differentially expressed genes.
Figure 5. DP cell differentiation trajectory and dynamic changes in the branch 1 gene. (a) Construction of the DP pseudotime differentiation trajectory; (b) changes in the expression of LEF1, MSX1, and HOXC13 in pseudotime; (c) heatmap displaying the branch 1 expression pattern; (d) GO enrichment analysis of differentially expressed genes.
Animals 13 03818 g005
Animals 13 03818 g006
Figure 6. Analysis of DP cell second branch cell fate. (a) Heatmap displaying the branch 2 expression pattern; (b) GO enrichment analysis of differentially expressed genes; (c) expression of the IGFBP3 and IGFBP7 genes in different stages of yak hair follicle cycle development; (d) changes in the expression of IGFBP3 and IGFBP7 in pseudotime.
Figure 6. Analysis of DP cell second branch cell fate. (a) Heatmap displaying the branch 2 expression pattern; (b) GO enrichment analysis of differentially expressed genes; (c) expression of the IGFBP3 and IGFBP7 genes in different stages of yak hair follicle cycle development; (d) changes in the expression of IGFBP3 and IGFBP7 in pseudotime.
Animals 13 03818 g006
Animals 13 03818 g007
Figure 7. Analysis of the cell fate of the third branch of DP cells.
Figure 7. Analysis of the cell fate of the third branch of DP cells.
Animals 13 03818 g007
Animals 13 03818 g008
Figure 8. Reasoning of intercellular communication in hair follicles. (a) Figures of the number and weight of cell–cell communication interactions of identified cell types: (b) intercellular communication between fibroblast types; (c) signaling pathway network diagram; (d) signaling-pathway-mediated intercellular communication; (e) output and input signal recognition of different cell groups. (f) distribution map of signal gene expression.
Figure 8. Reasoning of intercellular communication in hair follicles. (a) Figures of the number and weight of cell–cell communication interactions of identified cell types: (b) intercellular communication between fibroblast types; (c) signaling pathway network diagram; (d) signaling-pathway-mediated intercellular communication; (e) output and input signal recognition of different cell groups. (f) distribution map of signal gene expression.
Animals 13 03818 g008
Animals 13 03818 g009
Figure 9. Immunofluorescence analysis of key marker proteins in yak hair follicles.
Figure 9. Immunofluorescence analysis of key marker proteins in yak hair follicles.
Animals 13 03818 g009
Table 1. Single-cell data set quality statistics.
Table 1. Single-cell data set quality statistics.
SampleTelogenAnagen
Estimated number of cells14,57312,000
Median genes per cell1261927
Genome reads mapped87.6%93.3%
Valid barcodes97.7%97.1%
Total genes17,79618,420
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zheng, Q.; Ye, N.; Bao, P.; Wang, T.; Ma, C.; Chu, M.; Wu, X.; Kong, S.; Guo, X.; Liang, C.; et al. Interpretation of the Yak Skin Single-Cell Transcriptome Landscape. Animals 2023, 13, 3818. https://doi.org/10.3390/ani13243818

AMA Style

Zheng Q, Ye N, Bao P, Wang T, Ma C, Chu M, Wu X, Kong S, Guo X, Liang C, et al. Interpretation of the Yak Skin Single-Cell Transcriptome Landscape. Animals. 2023; 13(24):3818. https://doi.org/10.3390/ani13243818

Chicago/Turabian Style

Zheng, Qingbo, Na Ye, Pengjia Bao, Tong Wang, Chaofan Ma, Min Chu, Xiaoyun Wu, Siyuan Kong, Xian Guo, Chunnian Liang, and et al. 2023. "Interpretation of the Yak Skin Single-Cell Transcriptome Landscape" Animals 13, no. 24: 3818. https://doi.org/10.3390/ani13243818

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop