Identification of the Key Genes Associated with the Yak Hair Follicle Cycle

The development of hair follicles in yak shows significant seasonal cycles. In our previous research, transcriptome data including mRNAs and lncRNAs in five stages during the yak hair follicles (HFs) cycle were detected, but their regulation network and the hub genes in different periods are yet to be explored. This study aimed to screen and identify the hub genes during yak HFs cycle by constructing a mRNA-lncRNA co-expression network. A total of 5000 differently expressed mRNA (DEMs) and 729 differently expressed long noncoding RNA (DELs) were used to construct the co-expression network, based on weighted genes co-expression network analysis (WGCNA). Four temporally specific modules were considered to be significantly associated with the HFs cycle of yak. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the modules are enriched into Wnt, EMC-receptor interaction, PI3K-Akt, focal adhesion pathways, and so on. The hub genes, such as FER, ELMO1, PCOLCE, and HOXC13, were screened in different modules. Five hub genes (WNT5A, HOXC13, DLX3, FOXN1, and OVOL1) and part of key lncRNAs were identified for specific expression in skin tissue. Furthermore, immunofluorescence staining and Western blotting results showed that the expression location and abundance of DLX3 and OVOL1 are changed following the process of the HFs cycle, which further demonstrated that these two hub genes may play important roles in HFs development.


Introduction
The yak (Bos grunniens) is a key and symbolic species on the Tibetan Plateau. To resist the arctic-alpine and hypoxic environments, in addition to being covered with coarse wool in yak skin, a layer of undercoat (yak cashmere) grows on the bottom. Yak cashmere is a valuable textile material with fine fibers, similar to cashmere, and it has the characteristics of soft texture and warm performance, compared with coarse wool. This double-coated structure may be the result of yak adapting to the alpine environment for a long time. Coarse wool is derived from primary hair follicles (PHFs) and cashmere is derived from secondary hair follicles (SHFs). In yak, the scapular region, back, and side are prolific with cashmere, the abdomen is mainly coarse hair. The hair of yak is important for pastoralist's living materials and economic benefit in the Tibetan Plateau. Similar to cashmere goats, the SHFs of yak also undergo a clear seasonal circulation. A growth cycle mainly consists of growth (anagen), regression (catagen), and a rest stage (telogen) [1].
The hair follicles maintain its normal periodic activity is the result of response to multiple signaling molecules and their complex interactions, these signals may differently express in stages or specifically express in skin tissue. Over the past several decades, a large of signals involving in HFs development have been revealed and studied extensively.

Co-expression Network Construction and KEGG Enrichment Analysis
The R package WGCNA V1.64.1 was used to construct lncRNA-mRNA co-expression network based on the FPKM data of DEMs and DELs [21]. A signed weighted correlation network was built by calculating gene expression similarity to create an adjacency matrix. In this study, we selected a soft threshold power, β = 9, to establish the co-expression network, so that the correlation between genes conformed to the scale-free network distribution. One step network module was constructed by merging genes with highly similar coexpression patterns into modules. The dynamic tree cut algorithm-with parameters set as: minModeSize = 30, mergeCutHeight = 0.25-was employed to cut the hierarchal clustering tree. Gene modules were distinguished as different colors, and were indicated by different branches on the clustering tree. The eigengenes of each module was determined. To identify temporal specific modules, the correlations of these modules with the different developmental stages in the yak hair cycle were investigated using WGCNA package. The module eigengenes closely related to specific time points were screened as temporally specific modules. Protein-coding genes of each key module were submitted to DAVID for KEGG pathway enrichment analysis. Bos mutus was used as a background list in DAVID. The default parameters of the number of genes greater than 2 and p value < 0.1 were used as thresholds to output the KEGG pathway enrichment results of different modules. Top 20 KEGG pathways of different modules were visualized using ggplot2 R package.

Analysis of Temporally Specific Co-expression Network
The co-expression networks closely related to distinct stages were, respectively, analyzed and visualized by Cytoscape (version 3.5.1). Cytohubba was applied to calculate the topological parameters of the nodes and screen the key nodes identified as hub genes. Meanwhile, protein-protein interaction (PPI) network analyses were used to analyze the interaction of interested eigengenes into the STRING database (https://string-db.org, accessed on 6 June 2021). The medium confidence of the PPI network was set as the interaction score > 0.4 and visualized by Cytoscape.

RNA Extraction, Real-Time Quantitative PCR, and PCR Detection
Total RNA of the tissues was isolated using the Trizol reagent (Invitrogen, Wilmington, NC, USA). The RNA concentration and quality were evaluated using a NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The first-strand cDNA was synthesized using the PrimeScriptTM RT reagent kit (Takara, Dalian, China) according to the manufacturer's instructions. The RT-qPCR was performed using TB GreenTM Premix Ex Taq II (Takara, Dalian, China), and the following reaction conditions were used: 95 • C for 3 min, followed by 40 cycles of 95 • C for 10 s, and 60 • C for 30 s. The relative expression of genes was normalized using GAPDH and analyzed with the 2 −∆∆Ct method. PCR amplification was performed using the 2 × EasyTaq PCR SuperMix (TransGen Biotech, Beijing, China). The PCR procedure was as follows: 95 • C for 3 min, followed by 35 cycles of 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 20 s, and a final extension step was 72 • C for 5 min. β-actin and GAPDH were amplificated as two internal references. PCR amplification products were detected by 1.5% agarose gel. The primers used for RT-qPCR and PCR are listed in Table S3 and designed with oligo 6.

Western Blotting
Total protein of yak skin tissues was extracted with RIPA Lysis Buffer containing protease inhibitor and protein phosphatase (Beyotime, Shanghai, China), and protein concentrations were determined by BCA Protein Assay Kit (Beyotime, Shanghai, China). The protein expression levels of DLX3 and OVOL1 in different developmental stages were quantified using a capillary-based "WES" Simple Western System (ProteinSimple, San Jose, CA, USA), according to the manufacturer's instructions. The primary antibody information was as follows: anti-OVOL1 (PA5-41480; ThermoFisher; 1:50), anti-DLX3 (PA5-101065; ThermoFisher; 1:50), and GAPDH (10494-1-AP, Proteintech; 1:1000). The relative protein expression was calculated and analyzed based on the gel-like images produced by the Compass for SW software (Version 4.0, Protein Simple, SanJose, CA, USA).

Statistical Analysis
All the experiments were repeated three times. The statistical analyses were performed by GraphPad Prism 8.0 software and the data were presented as the mean ± SEM. Oneway ANOVA and the Dunnett's test were used for multiple comparisons. p < 0.05 was considered to be statistically significant.

Construction of mRNA-lncRNA Co-expression Network
Total 5000 differently expressed mRNAs (DEMs) (Table S1) and 728 differently expressed lncRNAs (DELs) ( Table S2) were identified and used to construct the mRNA-lncRNA co-expression network by WGCNA software package in R. Nine co-expression modules were established (Figure 1a), excluding the genes in grey module that cannot be assigned to co-expression module. Number of the genes in different modules varied widely, from 45 genes in the magenta module to 1599 in turquoise module (Table S4). An eigengene adjacency heatmap was generated to explore the correlations between the modules and the dendrogram branches with positively correlated eigengenes were grouped together ( Figure 1b). Then, the heatmap was established for the correlation analysis between the stages of HFs development and module eigengenes to identify the temporally specific modules. As shown in Figure 1c, each stage was related to a module that has a high correlation with it. The yellow module has a high correlation with Jan. (cor = 0.85, p = 4 × 10 −5 ); the turquoise module is closely associated with Mar. (cor = 0.74, p = 0.002); the red module is closely associated with Jun. (cor = 0.79, p = 5 × 10 −4 ); the blue module is significantly associated with Oct. (cor = 0.69, p = 0.004). Therefore, the four modules were screened as temporally specific modules for the further functional analysis. Due to Genes 2022, 13, 32 5 of 14 the black module closely associated with Aug., containing a smaller number of genes that failed in the functional enrichment analysis, Aug. was excluded. = 4 × 10 −5 ); the turquoise module is closely associated with Mar. (cor = 0.74, p = 0.002); the red module is closely associated with Jun. (cor = 0.79, p = 5 × 10 −4 ); the blue module is significantly associated with Oct. (cor = 0.69, p = 0.004). Therefore, the four modules were screened as temporally specific modules for the further functional analysis. Due to the black module closely associated with Aug., containing a smaller number of genes that failed in the functional enrichment analysis, Aug. was excluded.

Function Enrichment Analysis
To understand the biological functions of the temporally specific co-expression modules, KEGG enrichment analysis was performed by submitting the mRNAs in key modules to DAVID database (Tables S5-S8). KEGG enrichment results showed that the pathways, such as Basal cell carcinoma, Melanogenesis, Tight junction, and Wnt signaling pathway, were enriched in the blue module that related to Oct. (anagen), these pathways were reported to highly express anagen or involve induction of anagen hair follicle differentiation. The TNF signaling pathway, platelet activation, ECM-receptor interaction and PI3K-Akt signaling pathway, the pathways that involved in promoting hair follicle neogenesis, and those inducing folliculogenesis, were enriched in turquoise and red modules, the two modules are closely associated with Mar. and Jun., respectively. Mar. is in the telogen, and Jun. is in the transition of telogen to anagen during HFs cycle. In the yellow module, the cAMP signaling pathway, the MAPK signaling pathway, and the Ras signaling pathway were enriched, which were reported as involving in the induction of early catagen transition or hair loss. Furthermore, focal adhesion was enriched in several modules, including red, turquoise, and blue, implying that the pathway may play an essential role in hair follicles development ( Figure 2a). Meanwhile, the scatterplot of gene significance (GS) for each period vs. module membership (MM) of their associated modules were plotted (Figure 2b). The results showed that all the GS and MM exhibit very significant correlation, suggesting that hub genes in the key modules also tend to be highly correlated with their closely related stages, which also demonstrated the reliable biological significance of these co-expression modules.

Function Enrichment Analysis
To understand the biological functions of the temporally specific co-expression modules, KEGG enrichment analysis was performed by submitting the mRNAs in key modules to DAVID database (Tables S5-S8). KEGG enrichment results showed that the pathways, such as Basal cell carcinoma, Melanogenesis, Tight junction, and Wnt signaling pathway, were enriched in the blue module that related to Oct. (anagen), these pathways were reported to highly express anagen or involve induction of anagen hair follicle differentiation. The TNF signaling pathway, platelet activation, ECM-receptor interaction and PI3K-Akt signaling pathway, the pathways that involved in promoting hair follicle neogenesis, and those inducing folliculogenesis, were enriched in turquoise and red modules, the two modules are closely associated with Mar. and Jun., respectively. Mar. is in the telogen, and Jun. is in the transition of telogen to anagen during HFs cycle. In the yellow module, the cAMP signaling pathway, the MAPK signaling pathway, and the Ras signaling pathway were enriched, which were reported as involving in the induction of early catagen transition or hair loss. Furthermore, focal adhesion was enriched in several modules, including red, turquoise, and blue, implying that the pathway may play an essential role in hair follicles development ( Figure 2a). Meanwhile, the scatterplot of gene significance (GS) for each period vs. module membership (MM) of their associated modules were plotted ( Figure 2b). The results showed that all the GS and MM exhibit very significant correlation, suggesting that hub genes in the key modules also tend to be highly correlated with their closely related stages, which also demonstrated the reliable biological significance of these co-expression modules. Genes 2022, 13, x FOR PEER REVIEW 6 of 14

Screening of The Hub Genes in Temporally Specific Modules
To further analyze the potential regulatory relationship and screening the hub genes in the co-expression modules, the co-expression mRNA-lncRNA networks of the yellow, turquoise, red, and blues modules were analyzed, respectively. Genes of the top 100 degree of connectively in the different module networks were calculated and visualized using Cytohubba tab in Cytoscape. In the co-expression networks, the larger the node is, the higher connectivity in the network can be regarded as a potential hub regulator ( Figure  3). The results showed that FER, UGGT1, BDP1, and MED13 were identified to be the hub genes in the yellow module ( Figure 3a). ELMO1, CPEB1, and PID1 were located the hub position in the turquoise network (Figure 3b). KDELR3 and PCOLCE were the hub genes in the red module (Figure 3c). In the co-expression network of the blue module, multiple transcription factors, including HOXC13, FOXN1, MSX2, CUX1, DLX3, and OVOL1, and keratins, such as KRT35, KRT32, and KRT82, which were the crucial regulators or structure proteins for the HFs development, were revealed, and had a higher connectivity in the networks (Figure 3d). All the lncRNAs were at the relative lower degree in the networks. In general, this section screened and visualized the hub genes in different period of HFs development in yak.

Screening of The Hub Genes in Temporally Specific Modules
To further analyze the potential regulatory relationship and screening the hub genes in the co-expression modules, the co-expression mRNA-lncRNA networks of the yellow, turquoise, red, and blues modules were analyzed, respectively. Genes of the top 100 degree of connectively in the different module networks were calculated and visualized using Cytohubba tab in Cytoscape. In the co-expression networks, the larger the node is, the higher connectivity in the network can be regarded as a potential hub regulator ( Figure 3). The results showed that FER, UGGT1, BDP1, and MED13 were identified to be the hub genes in the yellow module ( Figure 3a). ELMO1, CPEB1, and PID1 were located the hub position in the turquoise network (Figure 3b). KDELR3 and PCOLCE were the hub genes in the red module (Figure 3c). In the co-expression network of the blue module, multiple transcription factors, including HOXC13, FOXN1, MSX2, CUX1, DLX3, and OVOL1, and keratins, such as KRT35, KRT32, and KRT82, which were the crucial regulators or structure proteins for the HFs development, were revealed, and had a higher connectivity in the networks (Figure 3d). All the lncRNAs were at the relative lower degree in the networks. In general, this section screened and visualized the hub genes in different period of HFs development in yak.  The networks are represented with their corresponding module colors, including the yellow, turquoise, red, and blue module networks (a-d). Genes in the red module network are less than 100 that were all exhibited. Circle represents mRNAs, diamond represents lncRNAs, the V shape in the blue network represents the keratins. The size of nodes represents the connectivity of genes in the co-expression network-the larger the nodes, the higher the connectivity and the more critical it is in the network and could be used as hub genes.

PPI Analysis and Expression Characteristics Detection of The Hub Genes in the Blue Module
Blue module was closely related to Oct. (Figure 1c), a period of full anagen during yak HFs cycle. In the blue module network, it was found that multiple hub genes were reported as being involved in HFs development. To further investigate the relationship between these hub regulators, PPI analysis was performed of the blue network. A subnetwork consisting of a series of transcription factors, such as HOXC13, MSX2, LEF1, DLX3, and FOXN1, and keratins including KRT32, KRT35, and KRT82, that were involved in HFs development, were presented in the PPI network (Figure 4a), which further indicates that these key genes may interact synergistically, promoting the hair growth. Then, the tissue specificity of part of the hub genes in the blue network and the PPI subnetwork were detected by PCR method. PCR result showed that, in addition to keratins, the expressions of WNT5A, HOXC13, DLX3, FOXN1, and OVOL1 were also relatively specific in their expressions in skin ( Figure 4b); these five genes also have been reported to be crucial regulators in HFs development. Then, the Pearson correlation coefficient between mRNAs and DELs were analyzed (Table S9), the DELs with a Pearson correlation coefficient (r) > 0.9 or -0.9 and p value < 0.05 with the five key regulators were visualized. There Figure 3. The co-expression network of the genes in top 100 connectivity in the key four modules. The networks are represented with their corresponding module colors, including the yellow, turquoise, red, and blue module networks (a-d). Genes in the red module network are less than 100 that were all exhibited. Circle represents mRNAs, diamond represents lncRNAs, the V shape in the blue network represents the keratins. The size of nodes represents the connectivity of genes in the co-expression network-the larger the nodes, the higher the connectivity and the more critical it is in the network and could be used as hub genes.

PPI Analysis and Expression Characteristics Detection of The Hub Genes in the Blue Module
Blue module was closely related to Oct. (Figure 1c), a period of full anagen during yak HFs cycle. In the blue module network, it was found that multiple hub genes were reported as being involved in HFs development. To further investigate the relationship between these hub regulators, PPI analysis was performed of the blue network. A subnetwork consisting of a series of transcription factors, such as HOXC13, MSX2, LEF1, DLX3, and FOXN1, and keratins including KRT32, KRT35, and KRT82, that were involved in HFs development, were presented in the PPI network (Figure 4a), which further indicates that these key genes may interact synergistically, promoting the hair growth. Then, the tissue specificity of part of the hub genes in the blue network and the PPI subnetwork were detected by PCR method. PCR result showed that, in addition to keratins, the expressions of WNT5A, HOXC13, DLX3, FOXN1, and OVOL1 were also relatively specific in their expressions in skin ( Figure 4b); these five genes also have been reported to be crucial regulators in HFs development. Then, the Pearson correlation coefficient between mRNAs and DELs were analyzed (Table S9), the DELs with a Pearson correlation coefficient (r) > 0.9 or −0.9 and p value < 0.05 with the five key regulators were visualized. There were 27 DELs that were found closely associated with the five regulators, among which, four DELs were negatively correlated, and others were positively correlated (Figure 4c). The expression of lncRNAs were reported to be spatially and temporally specific [25]. Six lncRNAs were randomly selected for tissue specific detection, and the result showed that part of the lncRNAs relatively specific expressed in skin (Figure 4b). Moreover, the expression changes of WNT5A, HOXC13, DLX3, FOXN1, and OVOL1 in different stages of HFs cycle were detected using RT-qPCR, which could also represent the expression trend of the genes in the blue module. The result showed that the five genes were highly expressed in anagen (Aug. and Oct.) (Figure 4d), consisting with the general expression trend of the blue module, analyzed by WGCNA. These results further identified and characterized the hub genes in blue, and analyzed the lncRNAs which interacted with key mRNAs. were 27 DELs that were found closely associated with the five regulators, among which, four DELs were negatively correlated, and others were positively correlated (Figure 4c). The expression of lncRNAs were reported to be spatially and temporally specific [25]. Six lncRNAs were randomly selected for tissue specific detection, and the result showed that part of the lncRNAs relatively specific expressed in skin (Figure 4b). Moreover, the expression changes of WNT5A, HOXC13, DLX3, FOXN1, and OVOL1 in different stages of HFs cycle were detected using RT-qPCR, which could also represent the expression trend of the genes in the blue module. The result showed that the five genes were highly expressed in anagen (Aug. and Oct.) (Figure 4d), consisting with the general expression trend of the blue module, analyzed by WGCNA. These results further identified and characterized the hub genes in blue, and analyzed the lncRNAs which interacted with key mRNAs.

Spatial-Temporal Expression Analysis of DLX3 and OVOL1 during Yak HFs Cycle
In order to further identify the hub genes that may play crucial roles in HFs development, DLX3 and OVOL1 were selected to analyze the spatiotemporal expression in protein level during the HFs cycle. DLX3 and OVOL1 were found located a core position in the

Spatial-Temporal Expression Analysis of DLX3 and OVOL1 during Yak HFs Cycle
In order to further identify the hub genes that may play crucial roles in HFs development, DLX3 and OVOL1 were selected to analyze the spatiotemporal expression in protein level during the HFs cycle. DLX3 and OVOL1 were found located a core position in the blue co-expression network and their roles in HFs development of yak have been poorly studied. Immunofluorescent staining results showed that the expression abundance and location of Genes 2022, 13, 32 9 of 14 DLX3 and OVOL1 in protein level changed during yak hair cycle. DLX3 expression was detected in the outer root sheath throughout the hair follicle cycle. Besides, by anagen stage DLX3 was specifically expressed in matrix and hair shaft, and was presented in epidermis chain of catagen. In telogen, DLX3 is highly expressed in the telogen bulge (Figure 5a). The result suggested that DLX3 may play a role in the whole process of the HFs cycle, a phenomenon which could also be reflected in Figure 5c. OVOL1 expression was detected in the hair matrix in anagen, and is mainly presented in the inner root sheath of the lower part of hair follicle in early anagen, full anagen, and later anagen. In catagen and telogen, the expression of OVOL1 is weaker, implicating that OVOL1 mainly regulates the growth of hair follicles in anagen (Figure 5b), and this expression distribution could also be observed in Figure 5c. Then, WB was used to determine the expression of DLX3 and OVOL1 in Jan. (catagen), Apr. (telogen), and Sep. (anagen), the result was basically consistent with the result of immunofluorescent. Both DLX3 and OVOL1 are highly expressed in Sep. (anagen). Compared with OVOL1, DLX3 has a higher expression abundance in each stage, and the expression of OVOL1 is lowest in Jan. (catagen) (Figure 5d).
Genes 2022, 13, x FOR PEER REVIEW 9 of 14 blue co-expression network and their roles in HFs development of yak have been poorly studied. Immunofluorescent staining results showed that the expression abundance and location of DLX3 and OVOL1 in protein level changed during yak hair cycle. DLX3 expression was detected in the outer root sheath throughout the hair follicle cycle. Besides, by anagen stage DLX3 was specifically expressed in matrix and hair shaft, and was presented in epidermis chain of catagen. In telogen, DLX3 is highly expressed in the telogen bulge (Figure 5a). The result suggested that DLX3 may play a role in the whole process of the HFs cycle, a phenomenon which could also be reflected in Figure 5c. OVOL1 expression was detected in the hair matrix in anagen, and is mainly presented in the inner root sheath of the lower part of hair follicle in early anagen, full anagen, and later anagen. In catagen and telogen, the expression of OVOL1 is weaker, implicating that OVOL1 mainly regulates the growth of hair follicles in anagen (Figure 5b), and this expression distribution could also be observed in Figure 5c.

Discussion
The hair cycle of mammals includes anagen, catagen, and telogen. The periodic transformation of the hair cycle is regulated by multiple signaling molecules and complex interactions between the epidermis and dermis. In animal husbandry, the normal hair cycle is of great significance to the economic benefit of fiber-producing animals, such as cashmere goat, sheep, and angora rabbit. The seasonal cycling of secondary HFs allows yak to have a thicker coat in cold winters and molt in warmer seasons, which are crucial for yak's resistance to the alpine and harsh environments. In the present study, the differently expressed mRNAs and lncRNAs of yak skin at different stages during HFs cycling (Tables S1 and S2) was used to construct the mRNA-lncRNA co-expression network by WGCNA. There were four temporally specific modules that were found closely related to different HFs development stages, respectively (Figure 1c). These modules could be applied to screen and identify the hub genes and important signaling pathways of different periods of HFs cycle. KEGG analysis result showed that the signaling pathways enriched in the modules are consistent with the biological process of that module closely related period ( Figure 2a). Wnt, basal cell carcinoma, melanogenesis, and tight junction were enriched in the blue module. Wnt has been considered as a key signaling pathway that drives HFs anagen induction and is required for hair follicle growth and regeneration [26,27]. Basal cell carcinoma, melanogenesis, and tight junction were reported to be highly expressed in anagen or as involving the induction of anagen hair follicle differentiation [28][29][30]. ECMreceptor interaction and PI3K-Akt signaling pathways were enriched in both the red and the turquoise modules. EMC plays an important role in the hair follicle stem cell nicheregulating, bulge cell behavior and hair follicle regeneration [31,32]. PI3K-Akt signaling is essential for de novo hair follicle regeneration [33]. Mar. is the telogen of yak HFs cycle; Jun. is the transitional period from telogen to anagen. During the stage of telogen to early anagen, the hair follicle stem cells transform from a quiescent to activated state, and this is a crucial stage of hair follicle neogenesis. Besides, platelet activation and TNF signaling pathways that involved in inducing folliculogenesis and promoting of hair follicle neogenesis were also enriched in turquoise and red modules [34,35]. In the yellow module, the cAMP signaling pathway [36], the MAPK signaling pathway [37], and the Ras signaling pathway [38] were reported to be involved in induction of early catagen transition or hair loss, the functions coincide with the biological process of HFs in Jan. (catagen). Thus, KEGG enrichment analysis revealed the biological functions of the temporally specific modules, as well indicated the validity of our module construction.
In this study, hub genes of the key modules were screened. Genes located at the core of the network could be focused on as key regulatory genes. Fer is a hub gene located at the core of the yellow module network (Figure 3a), which was reported to be involved in the growth and metabolism of various epithelial tumors and was deemed to an important positive regulator of melanoma metastasis. A lack of Fer kinase in melanocytes emerged an in impaired activity of Wnt/β-catenin signaling pathway [39]. In the turquoise module (Figure 3b), the hub gene ELMO1 was reported playing an unexpected role in the clearance of apoptotic cells and in the maintenance of homeostasis [40]; CPEB1 inhibits epidermmesenchymal transformation and regulates cell cycles [41]. These functions of the hub genes may contribute to the maintenance and dynamics in telogen. PCOLCE and KDELR3 are the prominent hub genes in red module associated with Jun. (Figure 3c). PCOLCE was reported to be an important factor in the transition from telogen to anagen [42], and the expression of KDELR3 is associated with the development of melanocytes and the progression of melanoma [43]. Furthermore, several hub genes found in the blue module (Figure 3d), including MSX2, FOXN1 [8], HOXC13 [44], DLX3 [12], OVOL1 [45], LEF1 [46], TCHH, and PADI3 [47,48], are well known to be crucial regulators of hair follicle development and are located in the core of the network. In addition, more keratins, such as KRT82, KRT35, KRT35, KRT28, and KRT71, were showed in the blue module. The blue module was related to anagen (Oct.); anagen is the longest stage of the hair follicle development cycle, accounting for more than half of the whole growth cycle. The anagen in yak hair cycle is from about Jun. to Dec. Moreover, researchers are also likely to focus on growing hair follicle due to the intact morphology, which may lead to more genes involved in regulation of hair follicle development to be found in the blue module. The intuitive display of hub genes in different stages of the hair cycle would provide reference information for further mining of regulators related to the hair cycle and reveal the mechanism of hair follicle development.
Multiple regulators and keratins in the blue co-expression network (Figure 3d) were also found to constitute a PPI network (Figure 4a), which further indicates that these genes may interact with each other in hair follicle development. Among these key genes, in addition to several keratins, the regulators WNT5A, HOXC13, DLX3, FOXN1, and OVOL1 were also found to be specifically expressed in skin tissue (Figure 4b). The regulators WNT5A, HOXC13, and FOXN1 have been widely reported to play important roles in hair follicle development [8,9,49].
DLX3 and OVOL1 were two hub genes in the blue module and were specifically expressed in skin. DLX3 is a transcription factor that regulates epidermal, neural, and osteogenic cellular differentiation [12]. OVOL1 is a putative transcription factor and is mainly involved in hair formation and spermatogenesis [50]. In the present research, DLX3 was shown to be expressed in the outer root sheath, throughout the hair follicle cycle, and expressed in matrix and hair shaft in anagen, and in the epidermis chain in catagen (Figure 5a). Detection of the expression position of DLX3 in our study is consistent with previously reported data [12]. Our results intuitively and globally demonstrated the expression of DLX3 in yak hair follicle. Changes of DLX3 in expression abundance and location during the HFs cycle further evidenced that DLX3 is a crucial regulator in hair morphogenesis, differentiation, and cycling programs. OVOL1 was detected, mainly presenting in the inner root sheath of the lower part of the hair follicle in the early anagen, mid-anagen, and later anagen, and was expressed in the hair follicle matrix in mid-anagen ( Figure 5b); the result is consistent with the confirmed expression of OVOL1 in the inner root sheath of human hair follicle [45,51]. The OVOL1-OVOL2 axis was reported to be crucial for the normal development of hair follicles in murine. OVOL1-deficient mice indicate ruffled hair coat and hair abnormalities [50]. However, the role of OVOL1 in the hair follicle cycle is rarely reported. Our study of OVOL1 expression changes during yak hair follicle cycle will provide clues to reveal the specific function of OVOL1 in HFs development.

Conclusions
In summary, a lncRNA-mRNA co-expression network was constructed to identify the key pathways and the hub genes closely related to the development of yak HFs cycle. A series of well-known biological processes or pathways associated with HFs development were enriched by analyzing the temporally specific co-expression gene modules. The crucial genes in different stages of yak HFs cycle were screened by co-expression network analysis, combined with different molecular experimental methods. Our findings systematically elucidated the biological processes and important regulators during the development of yak HFs cycle, which would contribute to better understanding of molecular mechanism in the development of hair follicle and provide a useful reference information for molecular breeding in yak hair traits.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes13010032/s1, Table S1: Differentially expressed mRNAs used for co-expression analysis, Table S2: Differentially expressed lncRNAs used for co-expression analysis, Table S3: Primers used in quantitative PCR analysis and PCR, Table S4: WGCNA module membership for each gene and  lncRNA, Tables S5-S8: Results of gene enrichment analysis for each module, including yellow, turquoise, red, and blue modules, respectively, Table S9: Pearson correlation analysis between mRNAs and DELs.