Proteomic analysis of fetal skin by iTRAQ reveals molecular signals underlying Inner Mongolia Cashmere goat hair follicle initiation

Abstract

Cashmere goats produce wool from primary hair follicles (PHF)and cashmere from secondary hair follicles(SHF). Numerous genetic factors have been identified to affect HF formation, but the mechanism that controls the distinct morphogenesis of PHFs and SHFs remains unclear. To better describe goat’s hair follicle proteome, we applied iTRAQ technology to study protein expression profiling from nine skin samples of cashmere goats in three fetal periods including 45-day-old embryos (E45), 55-day-old embryos (E55), 65-day-old embryos (E65). A total of 1172 proteins identified are expressed in fetal skin. The numbers of differentially expressed proteins in E45 vs. E55, E45 vs. E65, and E55 vs. E65 were 103, 144, and 118, respectively. Sixty-three proteins were constantly increased through the whole embryonic developmental stages of placode formation. Keratin 1, 4, 5, 10, 13, 15, 17, and 75 belonging to the hair keratin family (KRTs) were identified during keratinocyte proliferation. Some proteins involved in these pathways probably have a role in the HF and might be the mediators of HF initiation in Cashmere goats. Our systematic investigation identified some signaling factors that drive the large-scale cellular rearrangements necessary for HF initiation, which provide an insight into the fundamental mechanisms of cashmere growth.

KEYWORDS:

• ITRAQ
• Cashmere goat
• fetal skin
• hair follicle
• proteomics
• clustering

Introduction

Fiber fineness is very important in the luxury textile industry for manufacturing high-quality textiles. Cashmere is obtained from the undercoat hair of several domestic mammals, e.g. goats, alpaca. China is the world’s largest producer of cashmere fibers and Cashmere goats play an important role in the development of China’s animal husbandry. Inner Mongolia Cashmere goat has a double coat consisting of coarse wool and cashmere. The coarse wool comes from the primary hair follicles (PHFs) and cashmere is derived from the secondary hair follicles (SHFs). The value of cashmere fleece is mainly determined by the numbers and density of HFs (Allain and Renieri 2010). HF organogenesis largely takes place during fetal skin development. And reciprocal inductive interactions between the epithelium and mesenchyme have major effects on the development of hair follicles during embryonic day (Hardy 1992). A series of signals are involved in the hair follicle morphogenesis, including epithelial thickening, hair papilla formation, and deepening of hair follicle epithelium to the dermis. Specific ectodermal cell clusters are directed to enlarge into epithelial placodes by a signal coming from the mesenchyme. The placodes then send a signal to the underlying mesenchyme that instructs fibroblasts to condense into primitive dermal papillae. Subsequently, dermal condensates encourage epithelial cell growth and differentiation (Fuchs and Raghavan 2002; Millar 2002; Mill et al. 2003). The primary hair placodes are first evident in the lateral surface ectoderm at embryonic day (E) 18.5. (Mill et al. 2003) Our histological study has demonstrated that embryonic skin epidermis formed between E45 and E55 and two types of HF had not yet started at E45 in cashmere goats. The primary hair follicle primordium firstly appears on embryo mid-flank skin at E55. The original of secondary hair follicle appeared only on the side of neck at day 65 (E65) and spread through whole body by 75 days (Zhang et al. 2007).

The importance of proteomics in animal science has increased significantly in recent years. Since the success of the genome projects, the vast of genes/proteins has been associated with the desired traits for the domestic animals. Proteomics is can provide the insights into protein expression, protein–protein interactions, or post-translational modifications during the biological process. This technique will assist researchers to investigate the cellular functions in whole. Proteomic methods have been used in porcine and bovine, poultry species (Almeida et al. 2015). For hair follicle morphogenesis, proteomics analysis has been done in mouse (Wang et al. 2016), human dermal papilla cell (Zhang et al. 2016), which provided some understanding in molecules involved in hair follicle development mechanisms. However, for cashmere goat, most studies have been focused on wool quality, breed diversity in specific high sulfur protein group (Flanagan et al. 2002; Xiaohui 2012; Guishan 2016; Li et al. 2018), wool cuticle (Koehn et al. 2009, 2010), and wool fiber morphogenesis (Plowman et al. 2015).

The generation of HF is established early during fetal life. The quantity and quality of cashmere are associated mainly with the traits of skin follicle. It is necessary to understand skin follicle structure, development, and regulatory genes, to acquire goats with good quantity and quality cashmere via selective breeding. To the best of our knowledge, the effect of the proteome in skin HF development at fetal periods is rarely reported. Isobaric tag for relative and absolute quantitation (iTRAQ)-based technology has gained great popularity in quantitative proteomics applications because of its relatively high sensitivity and reproducibility. The iTRAQ technique is an effective method for the global comparison of protein expression levels in a relatively small number of samples. It is capable of simultaneously identifying and quantifying proteins from multiple samples. It allows the differential labeling of peptides from distinct proteomes and the multiplexing ability afforded by the iTRAQ reagent, which are available in four to eight different tags (Pierce et al. 2008). Therefore, we will explore differences in protein constitution among different developmental stages from Inner Mongolia Cashmere goat using the iTRAQ method. This study aims to investigate the expression profiles of skin proteins at three specific developmental stages: 45-, 55-, and 65-day-old embryos (E45, E55, and E65) using iTRAQ. Ultimately, we aim to identify the molecular basis of fetal HF morphogenesis in cashmere goats.

Materials and methods

Ethics statement

This study was approved by the Key Laboratory of Animal Genetics, Breeding and Reproduction, Inner Mongolia Agricultural University, China. All surgery was performed under pentobarbital anesthesia and all efforts were made to minimize suffering. The procedure was approved by the Special Committee on Scientific Research and Academic Ethics of Inner Mongolia Agricultural University, responsible for the approval of biomedical research ethics of Inner Mongolia Agricultural University (Approval No. [2020] 056). No specific permissions were required for these activities, and no endangered or protected species were involved.

Tissue samples

The goats were housed at the Aerbasi White Cashmere Goat Breeding Farm (Inner Mongolia, China). Multiparous ewes aged four to five years were grazed under the same environment. Feeding and water were free access. In total, nine goat embryos were acquired by cesarean surgeries in the 45-, 55-, and 65-day-old pregnant goats with three biological repeats at each period. The skin samples were collected from the lateral body of each fetal goat and immediately frozen in liquid nitrogen for storage and transport until used for protein isolation.

Protein digestion, iTRAQ labeling, and strong cationex change

The protein digestion as described by Magdalena Luczak et al. (2016).

iTRAQ labeling was conducted following the manufacturer’s recommendations (Applied Biosystems, Foster City, CA, USA). Each group of samples was labeled as (45 days embryo)-114, (55 days embryo)-115, (65 days embryo)-116, and (the reference)-117, and were multiplexed and vacuum dried. The reference was composed of an equal amount of the different labeled samples.

iTRAQ-labeled peptides were fractionated on SCX cartridge system (AB Sciex, USA). The peptides were eluted with a sequence of buffer at a flow rate of 1 ml/min: (1) a gradient of 0%–10% buffer B (500 mM KCl, 10 mM KH2PO4 in 25% of ACN, pH 2.7) for 2 min, (2) 10–20% buffer B for 25 min, (3) 20%–45% buffer B for 5 min, and (4) 50%–100% buffer B for 5 min. Monitoring of the elution process was performed by measuring the absorbances at a wavelength of 214 nm. And fractions were collected with a rate of 1 min. Finally, all the collected fractions were pooled and divided into 10 parts which were desalted on C18 Cartridges (Empore™ SPE Cartridges C18 (standard density), bed I.D. 7 mm, volume 3 ml, Sigma). Every fraction was vacuum centrifuged for concentration and reconstituted in 40 µl of 0.1% (v/v) trifluoroacetic acid. All samples were stored at −80°C until LC-MS/MS analysis.

LC-MS/MS analysis

LC-MS/MS Analysis based on the detailed description by Magdalena Luczak et al. (2016).

Protein identification and quantification

MS/MS spectra were searched using MASCOT engine (Matrix Science, London, UK; version 2.2) against the UniprotKB Caprinae database (downloaded at April 2014; 30,424 sequences). For protein identification, the following options were used: Peptide mass tolerance = 20 ppm, MS/MS tolerance = 0.1 Da, Enzyme = Trypsin, Missed cleavage = 2, Fixed modification: Carbamidomethyl (C), iTRAQ 4plex (K), iTRAQ 4plex (N-term), Variable modification, Oxidation(M), and Decoy database pattern = Reverse.

The MASCOT search results for each SCX elution were further processed using the Proteome Discoverer1.4 (Thermo) (http://www.proteomics.ac.cn/). Assembling proteins were identified by Build Summary program based on a target-decoy search in shotgun proteomics. The criterion of reported data was set on 99% confidence for protein identification, which was determined by false discovery rate (FDR) ≤ 1%. FDR = N(decoy)*2/((N(decoy) + N(target)).

The ratios of proteins were determined by the programs Isobaric Labeling Multiple File Distiller and Identified Protein iTRAQ Statistic Builder; the mixture was used as reference. The expressions of proteins were calculated by the ratios of proteins normalized by a reference.

The protein–protein interaction (PPI) information of the studied proteins was retrieved from IntAct molecular interaction database (http://www.ebi.ac.uk/intact/) by their gene symbols or STRING software (http://string-db.org/). The results were downloaded in the XGMML format and imported into Cytoscape software (http://www.cytoscape.org/, version 3.2.1) to visualize and further analyze functional protein–protein interaction networks. Furthermore, the degree of each protein was calculated to evaluate the importance of the protein in the PPI network. (Shannon et al. 2003; Ono et al. 2014)

Bioinformatics annotation

Differentially expressed proteins (DEPs) were subjected to KEGG pathway and GO annotation using cluego in Cytoscape 3.7.2 (http://apps.cytoscape.org/ apps/cluego). The studied proteins were blasted against the KEGG database (http://www.KEGG.jp) and were subsequently mapped to pathways in KEGG. GO functional enrichment analysis was carried out using the GO database (http://geneontology.org/).

Figure 1. Morphological observation of hair follicles across three developmental stages in cashmere goat. (A–C) Skin at embryonic day 45 (E45), 55 (E55), 65 (E65), respectively. (A) Forty-five days of embryonic period: dermal signal initiation; The induction of hair follicles is initiated around E45 and characterized by crowded epidermal keratinocytes, hair placode, and dermal condensate. (B) Fifty-five days of embryonic period: primary hair follicle genesis; the number of primary follicles increases significantly. (C) Sixty-five days of embryonic period: secondary hair follicle genesis, secondary follicles (SFs) start to form.

Results

Phenotype of three stages of embryo

The phenotypes of the three stages of the embryo are shown in Figure 1.

Identified proteins analysis

iTRAQ labeling technology coupled with LC-MS/MS and Proteome Discoverer1.4 (thermo) was employed for protein identification and quantification in the three groups of 4-plex iTRAQ experiments. In total, 7753 unique peptides that corresponding to 1172 proteins were identified in the skin samples from three periods of cashmere goat fetuses (Supplementary Table S1). Of the proteins identified, 94.0% were represented by one to five peptides (Figure 2(A)), and the protein molecular weight was mainly distributed between 10 and 110 kDa (Figure 2(B)). The majority of the proteins were identified with high peptide coverage (Figure 2(C)) and the fold changes of the 55/45 were mostly close to 1 (Figure 2(D))

Figure 2. (A) Distribution of proteins containing different numbers of identified peptides, (B) distribution of the identified proteins among the different molecular weight classes, (C) coverage of proteins by the identified peptides, (D) distribution of protein abundance ratio.

Differentially expressed proteins in E45, E55, and E65

To better examine the biological mechanism of HF morphogenesis, it is important to identify the differentially expressed proteins among each developmental stage. Each developmental stage contained three biological replicates, and every developmental stage was considered as a group. Comparing the three paired groups (E55 vs. E45, E65 vs. E45, and E65 vs. E55), we detected eight overlapped proteins from the three groups. The numbers of the differentially expressed proteins in E55 vs. E45, E65 vs. E45, and E65 vs. E55 were 103, 144, and 118, respectively, with filtering thresholds of Fold change >1.2 or Fold change <0.83 and P-value ≤.05. Furthermore, E65 vs. E45 had the most differentially expressed proteins (Number = 144). Forty-seven common differentially expressed proteins between E55 vs. E45 and E65 vs. E45 were found. Aside from the 47 common differential proteins, 97 differentially expressed proteins were identified from E45 vs. E65, suggesting that more proteins were differentially regulated between E55 and E65 during embryonic placode formation.

Expression patterns and hierarchical cluster analysis during the three developmental stages

There are 188 expressed proteins significantly different (P < .05) in three hair follicle developmental stages (Table S2); 159 proteins were obtained according to the fold change greater than 1.2 or less than 0.83. To better examine the biological mechanism of HF morphogenesis, a heat map of 159 proteins illustrates distinct protein profiles by HF developmental stage and rationality of sampling (Figure 3). Moreover, 63 proteins were constantly increased through the whole embryonic developmental stages of placode formation. Among these proteins, we identified Keratin 1, 4, 5, 10, 13, 15, 17, and 75, which belong to the hair keratin family (KRTs) during keratinocyte proliferation in cashmere goats. These were enriched in the cuticle and may function in HF morphogenesis. The following proteins were also constantly increased extracellular matrix (ECM) components, TGFB1, and Decorin (a part of connective tissue) (Table 1). Forty-five proteins were constantly decreased, and 41 proteins showed a decreasing trend before they rose in the three developmental stages. Our results indicate that these differentially expressed proteins may exert fundamental function in goat embryonic HF initiation. A total of eight proteins were identified as differentially expressed proteins in all three groups, which may contribute to HF initiation and development (Figure 4). The information of these eight proteins has shown in Table 2 in the manuscript.

Figure 3. Hierarchical cluster analysis of protein expression profiles from nine skin samples with 159 proteins.

Table 1. Selected proteins constantly increased through the whole embryonic developmental stages of placode formation in cashmere goat skin.

Accession	Symbol	Description
Hair keratin family
O77727	KRT15	Keratin, type I cytoskeletal 15 OS = Ovis aries GN = KRT15 PE = 2 SV = 1 – [K1C15_SHEEP]
W5Q6E8	KRT75	Uncharacterized protein OS = Ovis aries GN = KRT75 PE = 3 SV = 1 – [W5Q6E8_SHEEP]
W5Q665	KRT13	Uncharacterized protein OS = Ovis aries GN = KRT13 PE = 3 SV = 1 – [W5Q665_SHEEP]
W5Q5U5	KRT4	Uncharacterized protein OS = Ovis aries GN = KRT4 PE = 3 SV = 1 – [W5Q5U5_SHEEP]
W5Q611	KRT1	Uncharacterized protein (Fragment) OS = Ovis aries GN = KRT1 PE = 3 SV = 1 – [W5Q611_SHEEP]
W5Q687	KRT5	Uncharacterized protein OS = Ovis aries GN = KRT5 PE = 3 SV = 1 – [W5Q687_SHEEP]
W5Q160	KRT10	Uncharacterized protein OS = Ovis aries GN = KRT10 PE = 3 SV = 1 – [W5Q160_SHEEP]
U3MXA9	KRT17	Keratin 17 OS = Ovis aries GN = KRT17 PE = 2 SV = 1 – [U3MXA9_SHEEP]
Annexins
W5PSZ5	ANXA1	Annexin OS = Ovis aries GN = ANXA1 PE = 3 SV = 1 – [W5PSZ5_SHEEP]
A2SW69	ANXA2	Annexin A2 OS = Ovis aries GN = ANXA2 PE = 1 SV = 1 – [ANXA2_SHEEP]
W5PEC3	ANXA6	Annexin OS = Ovis aries GN = ANXA6 PE = 3 SV = 1 – [W5PEC3_SHEEP]
Extracellular matrix (ECM) components
W5Q5H8	FGA	Fibrinogen alpha chain OS = Ovis aries GN = FGA PE = 4 SV = 1 – [W5Q5H8_SHEEP]
W5Q4M3	COL5A2	Uncharacterized protein OS = Ovis aries GN = COL5A2 PE = 4 SV = 1 – [W5Q4M3_SHEEP]
Q9MZW2	COL6A1	Type VI collagen alpha 3 chain (Fragment) OS = Ovis aries GN = COL6A1 PE = 2 SV = 1 – [Q9MZW2_SHEEP]
W5PTN6	COL6A2	Uncharacterized protein OS = Ovis aries GN = COL6A2 PE = 4 SV = 1 – [W5PTN6_SHEEP]
W5QCP9	COL6A3	Uncharacterized protein OS = Ovis aries GN = COL6A3 PE = 4 SV = 1 – [W5QCP9_SHEEP]
W5Q4S0	COL3A1	Uncharacterized protein OS = Ovis aries GN = COL3A1 PE = 4 SV = 1 – [W5Q4S0_SHEEP]
Growth factor and associated molecules
W5Q4S0	COL3A1	Uncharacterized protein OS = Ovis aries GN = COL3A1 PE = 4 SV = 1 – [W5Q4S0_SHEEP]
W5Q0F3	TGFBI	Uncharacterized protein (Fragment) OS = Ovis aries GN = TGFBI PE = 4 SV = 1 – [W5Q0F3_SHEEP]
Q9TTE2	DCN	Decorin OS = Ovis aries GN = DCN PE = 2 SV = 1 – [PGS2_SHEEP]

Table 2. Specific information of these eight proteins in the form of tables.

	Protein expression
Description	45–1	45–2	45–3	55–1	55–2	55–3	65–1	65–2	65–3	ANOVA p-value
KRT15	0.8453	0.6505	0.7181	1.3406	0.9549	0.8693	1.2997	1.2695	1.0076	.0235
KRT75	0.9311	0.9189	0.8163	0.9755	0.6965	0.8107	1.7067	1.4699	1.5065	.0019
KRT13	0.2857	0.4753	0.4084	0.9288	0.6141	0.9441	1.6172	2.2305	2.0617	.0053
KRT4	0.5475	0.6275	0.6620	0.9492	0.7447	0.9118	1.4046	1.9674	1.7361	.0093
KRT1	0.7870	0.8393	0.7176	0.8253	0.8229	0.9061	2.0776	1.9124	1.6516	.0026
KRT5	0.5297	0.5369	0.5632	0.8161	0.5326	0.8755	2.4429	2.1586	1.9363	.0019
KRT10	0.7557	0.6805	0.6544	0.7807	0.7984	0.7595	1.7205	1.8384	1.5921	.0003
KRT17	0.6063	0.7183	0.6827	0.8900	0.7321	0.8426	1.6998	1.9234	1.8364	.0003

Functional distribution of differentially expressed proteins

Among 188 differentially expressed proteins (P < .05) in the three different hair follicle developmental stages, 187 of these were assigned with one or more GO terms for biological process, molecular function, and cellular component according to gene ontology annotation (Figure 5). These proteins may have a more significant influence on the cell cycle and signal transduction. Among subcategories during HF morphogenesis, cell and cell part were prevalent in the main category of cellular components. Cellular process and binding were dominant in the main categories of biological process and molecular function.

Figure 4. Venn diagram depicting the overlap between proteins differentially expressed between E55 vs. E45, E65 vs. E55, and E65 vs. E45.

Figure 5. GO Analysis of differentially expressed Proteins in three developmental stages of hair follicles (P < .05).

We further analyzed the188 differentially expressed proteins by KEGG signal pathways, which could directly indicate if the enriched proteins were related to the occurrence and development of skin follicles. A total of 97 proteins were assigned with one or more KEGG annotations and were assigned to 99 KEGG pathways. 11.34% of these proteins were classified into metabolic pathways, including carbon fixation, nitrogen metabolism, lipid metabolism, and contain the most annotated proteins. Seven homologous proteins (COL6A1, COL6A2, COL6A3, FN1, HSPG2, THBS1, TNC) in ECM–receptor interaction, 7 (CDC42, COL6A1, COL6A3, FLNB, FN1, THBS1, TNC) in Focal adhesion, 6 (COL6A1, COL6A3, FN1, HSP90AB1, THBS1, TNC) in PI3K-Akt signaling, 3 (CDC42, FLNB, NFATC1) in MAPK signaling, and 2 (DCN, THBS1) in TGF-beta were detected in HF morphogenesis. VEGF and p53 signaling pathways were also obtained based on the KEGG pathway analysis (Table 3).

Table 3. KEGG pathway analysis.

Pathway	Count	Proteins
Metabolic pathways	11	ALDH9A1, ANPEP, BCAT1, CKB, CKMT2, CS, GFPT1, ISYNA1, PHGDH, PTGES3, PTGIS
ECM–receptor interaction	7	COL6A1, COL6A2, COL6A3, FN1, HSPG2, THBS1, TNC
Focal adhesion	7	CDC42, COL6A1, COL6A3, FLNB, FN1, THBS1, TNC
PI3K-Akt signaling pathway	6	COL6A1, COL6A3, FN1, HSP90AB1, THBS1, TNC
MAPK signaling pathway	3	CDC42, FLNB, NFATC1
TGF-beta signaling pathway	2	DCN, THBS1
VEGF signaling pathway	1	CDC42
Wnt signaling pathway	1	NFATC1
p53 signaling pathway	1	THBS1

PPI network analysis

The PPI network for differentially expressed proteins had 128 nodes and 1068 interactions (Figure 6). The degrees of the nodes in the PPI network ranged between 1 and 71. Especially, HSP90AB1 (degree = 71), GAPDH (degree = 38), ANXA2 (degree = 19), ATP5B (degree = 12), ATP5A1 (degree = 11), and ATP5D (degree = 6) were with higher degrees in the PPI network (Table S3, Supporting Information).

Figure 6. (A) The protein–protein interaction network for the differentially expressed proteins. (B) The differentially expressed proteins are highlighted in yellow.

Discussion

The generation of HF is established early during fetal life. In the present study, we performed an iTRAQ-based proteomics to analyze the expression profiles of skin proteins at three specific developmental stages. iTRAQ allows rapid quantification of proteins of very low abundance in the skin at fetal periods. For this reason, we have systematically characterized the proteome profiles involved in the regulation of HF morphogenesis in the fetal period.

Three specific developmental stages are typical (E45, E55, and E65) during the time interval spanning initiation and development of HFs (Zhang et al. 2007). Our previous studies have shown the hair follicle morphogenesis process and development rule of cashmere goat in three embryonic periods. The original of primary hair follicle could be seen firstly over all goat skin during 55–65 day of fetal stage, the original of secondary hair follicle appeared only in the side of neck in day 65 and spread through whole body in 75 day. Skin hair follicle morphogenesis is not initiated embryo on day 45, followed by primary hair follicle begins to occur on day 55; lastly, secondary hair follicle genesis takes place on day 65 (Zhang et al. 2007).

This study using the iTRAQ approach provides the largest proteome profiles in embryonic cashmere goat to date, including 1172 proteins from three periods of cashmere goat fetuses. Many of the newly identified differentially expressed proteins belong to KRTs and KRTAPs, the major structural proteins of the hair fiber and sheath (Langbein and Schweizer 2005; Giesen et al. 2011). The occurrence and development of HFs are regulated by interactions between specific dermal fibroblasts and epithelial cells. Dermal fibroblasts are one of the main components of the skin and KRTs and KRTAPs are important structural proteins of skin keratinization. Epithelial cells and dermal fibroblasts maintain the normal structure and functions of HFs during maturation.

Molecular communication between the epidermis and the underlying mesenchyme, including changes in the expression/activity of a constantly growing number of cytokines, hormones, enzymes, neurotransmitters, and transcription factors, is a key prerequisite for HF morphogenesis. (Schmidt-Ullrich and Paus 2005; Mikkola 2007) These transformations are controlled by changes in the local signaling milieu (Geng et al. 2013), based on Growth factors, TGF-β1, Decorin, and COL3A1 were up-regulated expressed through the whole embryonic developmental stages of placode formation. Annexins are a family of calcium-dependent phospholipid-binding proteins. They are involved in the structural organization of the cell, cell differentiation, proliferation, apoptosis, and ion fluxes by interacting with various cell-membrane components. ANXA1 plays these roles by binding to other proteins such as epithelial growth factor receptors, actins, and integrin A4 (Moss and Morgan 2004). The most compelling molecules discovered to impact HF morphogenesis are ECM components. The extracellular matrix plays an important role in tissue and organ morphogenesis. The interactions of extracellular matrix with cells are mediated by transmembrane molecules such as integrins, proteoglycans, CD36, and other cell surface-associated components (Hubmacher and Apte 2013). ECM components, such as FGA, COL5A2, COL6A1, COL6A2, COL6A3, and COL3A1, are up-regulated expressed in the present study. The proteins associated with the focal adhesion, another important cell-to-matrix communication pathway, were also be found in the KEGG pathway analysis. Given the increased expression of keratins as well as the up-regulation of proteins involved in cell-extracellular matrix communication, it’s suggesting that goat cashmere HF morphogenesis involves increased keratin, cytoskeleton, and extracellular matrix creation process. However, many of the constantly increased proteins identified herein do not have a direct effect on HF morphogenesis. Further analyses will be necessary to understand how these proteins are involved in HF morphogenesis.

In our KEGG pathway-based analysis, several signaling pathways were found to be involved in the hair follicle development, including the ECM–receptor interaction, PI3K-Akt, MAPK, TGF-beta, VEGF, WNT, and p53 pathways (Table 4). This result is consistent with the previous study based on combined analysis of mRNA and miRNA data (Geng et al. 2013). Other species also suggested that multiple signaling pathways, together with specific transcription factors, are involved in hair follicle initiation and development (Kaplan and Holbrook 1994; Millar 2002; Mill et al. 2003; Schmidt-Ullrich and Paus 2005; Fuchs 2007; Kanehisa et al. 2008; Schneider et al. 2009; Plasari et al. 2010; Xu et al. 2013; Zhu et al. 2013; Jiang et al. 2014; Han et al. 2020). Therefore, it suggested that these pathways related to proteins probably have a role in the HF and might be the mediators of HF initiation in Cashmere goats.

Table 4. KEGG pathway analysis.

Pathway	Count	Proteins
Metabolic pathways	11	ALDH9A1, ANPEP, BCAT1, CKB, CKMT2, CS, GFPT1, ISYNA1, PHGDH, PTGES3, PTGIS
ECM–receptor interaction	7	COL6A1, COL6A2, COL6A3, FN1, HSPG2, THBS1, TNC
Focal adhesion	7	CDC42, COL6A1, COL6A3, FLNB, FN1, THBS1, TNC
PI3K-Akt signaling pathway	6	COL6A1, COL6A3, FN1, HSP90AB1, THBS1, TNC
MAPK signaling pathway	3	CDC42, FLNB, NFATC1
TGF-beta signaling pathway	2	DCN, THBS1
VEGF signaling pathway	1	CDC42
Wnt signaling pathway	1	NFATC1
p53 signaling pathway	1	THBS1

PPI network analysis revealed a close relationship among the differentially expressed proteins. Adenosine is a nucleoside that is naturally present in cells and is known to have various physiological functions. It has been previously reported that adenosine which acts as a potent regulator of hair growth, stimulates the growth of hair follicles by triggering the expression of some growth factors (Oura et al. 2008). However, the relationship between adenosine and growth factors is not fully understood. Our investigation provides insights into understanding the protein–protein interactions at work during cashmere hair follicle initiation and skin development.

Conclusions

Hair follicle morphogenesis requires a tightly coordinated balance of activities involving ECM–receptor interaction, PI3K-Akt, MAPK TGF-beta, WNT, and other signaling pathways. In the skin proteome analysis using iTRAQ, we identified 188 differentially expressed proteins and some signaling factors. They drive the large-scale cellular rearrangements necessary for cashmere HF morphogenesis. The 188 differentially expressed proteins include Keratin 1, 4, 5, 10, 13, 15, 17, 75, ANXA1, TGFB1, and Decorin are the key proteins found in this study and we will continue to pay attention to them in the future. These key proteins may be considered as potential candidate proteins for studying the dynamic epithelial–mesenchymal interactions at work during cashmere hair follicle initiation and skin development.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository (Ma et al. 2019) with the dataset identifier PXD035366.

Additional information

Funding

The reported work was supported by the National Natural Science Foundation of China (grant number 31560619) and the Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region (NMGIRT2322).