A Screen for Key Genes and Pathways Involved in High-Quality Brush Hair in the Yangtze River Delta White Goat

The Yangtze River Delta White Goat is the only goat breed that produces high-quality brush hair, or type III hair, which is specialized for use in top-grade writing brushes. There has been little research, especially molecular research, on the traits that result in high-quality brush hair in the Yangtze River Delta White Goat. To explore the molecular mechanisms of the formation of high-quality brush hair, High-throughput RNA-Seq technology was used to compare skin samples from Yangtze River Delta White Goats that produce high-quality hair and non high-quality hair for identification of the important genes and related pathways that might influence the hair quality traits. The results showed that 295 genes were expressed differentially between the goats with higher and lower hair quality, respectively. Of those genes, 132 were up-regulated, 62 were down-regulated, and 101 were expressed exclusively in the goats with high-quality brush hair. Gene Ontology and Metabolic Pathway Significant Enrichment analyses of the differentially expressed genes indicated that the MAP3K1, DUSP1, DUSP6 and the MAPK signaling pathway might play important roles in the traits important for high-quality brush hair.


Introduction
The Yangtze River Delta White Goat is the unique goat species that can produce high-quality meat, skin, and Type III hair. Brush hair is usually divided into three grades, Type I low-quality hair, Type II mid-quality hair, and Type III high-quality hair [1], which is specialized for making top-grade writing brushes and has many characteristics such as white color, straight peak, fine luster, and rich elasticity. The Yangtze River Delta White Goat, also called the Haimen goat, is found mainly in the Yangtze River delta plain, the central area of which is located in Haimen County, Nantong, Jiangsu Province [2]. It also has a laudatory name, brush hair goat, for its special production of hair. With a price up to 2,500~3,000 US dollar per kg, Type III hair is popular in Japan, Korea, Singapore, and other Southeast Asian countries influenced by Chinese culture and usually a traditional Chinese export in tight supply [3]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Type I hair produces coarse hair, its surface scales resemble fish scales, and its hair-shaft center contains a large number of medulla. Type II hair also produces coarse wool, its medullary cavity is small, with relatively few medulla. Type III hair is made up of cuticula and cortex. The hair-shaft center contains a small amount of medulla. Type III hair has the characteristics of fine wool, heterotypical hair, and coarse hair. The cuticula consist of circular scales, hybrid scales, and fish-scale scales. The upper part of the medullary cavity consists of transparent or spotty marrow. The middle of the medullary cavity is fine and contains fewer medulla. The under part of the medullary cavity has a large number of medulla [4]. Type III hair grows only on the back and neck ridge of the animal, whereas Type I and Type II hair are widely distributed [4] . Non high-quality brush hair refers to Type I and Type II hair. High-quality brush hair refers particularly to Type III hair, which is a unique product of the Yangtze River Delta White Goat [5].
The hair root from outside to inside is made up of the dermal sheath, the outer root sheath, and the inner root sheath (IRS). The IRS determines the character of the hair shaft. The IRS proliferates upward inside of the outer root sheath to support the growth of the fiber. In the sebaceous-glands layer, the IRS cells degenerate into pieces and finally disappear, releasing the hair fiber [6]. The production of brush hair is seasonal, which might be related to the growth cycle of the hair follicles [7]. The hair follicle growth cycle includes anagen, catagen, and telogen [8]. Cytokines play an important role in the regulation of the hair follicle growth cycle. Hair follicle stem cells have a certain life cycle, and the hair follicle subsequently enters catagen, the split termination [9]. The hair follicle growth and development process involves many genes and signaling molecules. Most of those signaling molecules belong to the Winglessrelated (Wnt) pathway, Sonic Hedgehog (Shh) pathway, MAPK pathway, and so on [10][11]. Previous studies have indicated that the MAPK pathway plays an important role in the growth and development of hair follicles [12][13][14]. The MAPK pathway can regulate cell growth, differentiation, and apoptosis and can participate in the regulation of hair-follicle growth and development [15][16]. In this study we used RNA-Seq technology to explore the important genes and pathways related to the formation of type III hair. The annotation of Unigenes from both the high hair quality and the non high hair quality phenotypes was obtained, and a gene pool containing all of the differentially expressed genes of the metabolic pathways was established to provide supporting evidence for the further study of the molecular mechanisms underlying the Type III hair traits in the Yangtze River Delta White Goat.

Ethics statement
The animal experiments were carried under the regulation of Yangzhou University Ethics Committee, and we have attached the Certificate on Animal experiments in this study (the certificate in the form of Supporting Information attached at the end of article), which is approved by Yangzhou University Ethics Committee. I confirm that the field studies did not involve endangered or protected species. This study location in the Haimen State Haimen Goat Farm and key laboratory for animal genetics and molecular breeding of Jiangsu Province, college of animal science and technology, Yangzhou University, and we have got their permission. I confirm that I received permission in advance of the study being conducted from your local ethics committee, and this permission for this specific study. The animal care and use were fully compliant with national standards of China. We followed the protocols and no animals were sacrificed for the purposes of this study. We followed the management and regulations for the goats. Rams were fed alone, and ewes fed together, and each goat occupied the living space of more than 1 square meters. The goat shed have the good condition of ventilation and penetrating light. We fed it once a day with free drinking water. We cleared the feces and dirt every day and disinfected the shed once a week. After operations, the skin cuts of the experimental animals were disinfected and cared well to restore.

Preparation of test samples
The experiments of goats came from the Haimen State Haimen Goat Farm. Six half-sib rams of 6~8 months were selected, including three high-quality hair and three non high-quality hair. From each goat, about 2 cm 2 of skin along the cervical spine was sampled and cryopreserved. The goats were given local anesthesia when the samples were taken. The six samples were divided into two groups: group A, comprising the samples from the three high hair quality goats, and group B, comprising the samples from the three non high hair quality goats.

Establishment and verification of the cDNA Library
First, the total RNA was extracted from the skin samples and tested for integrity. Then, a cDNA library was established and verified. The Illumina HiSeq TM 2000 system was used to sequence the transcript library. The three skin samples from the goats with the high hair quality traits were mixed for sample A, and the three skin samples from the goats with the non high hair quality traits were mixed for sample B. The samples were sequenced by computer, and the original sequencing data was stored in the FASTQ format. The sequencing work was performed by Shanghai Paisennuo Biotechnology Co. LTD. After the quality analysis of the original data, the low-quality sequences were removed, and the remaining high-quality sequences were analyzed further.

Data processing
The filtered sequences were spliced into transcripts from scratch (de novo assembly) using the Trinity software (Broad Institute). Then, each transcript were compared to the sequences in NR Protein, SwissProt, Kyoto Encyclopedia of Genes and Genomes (KEGG), and eggNOG databases using BLAST and to the Gene Ontology (GO) and Unigene annotations using Blas-t2GO (https://www.blast2go.com/). The KEGG annotation was used to obtain the pathway annotations of Unigene. Next, the Unigenes were analyzed using eggNOG.
Unigene expression and analysis of differentially expressed genes. Unigene expression was calculated by the Reads Per Kilobase per Million reads (RPKM) method [17] using the following formula: RPKM = 10 9 /NL. Where RPKM is the gene expression, C is the number of reads mapped to the gene, N is the total number of mapped reads, and L is the number of bases in the gene.
GO functional categories enriched with differentially expressed genes. GO is the international standard for the classification of gene functions. GO annotation includes entries in three main categories: Molecular Function, Cellular Component, and Biological Process. The GO annotations were obtained for each differentially expressed gene. The annotations were then analyzed using the WEGO software [18] to get a p-value for each GO entry. If the p-value was less than 0.5, then the GO entry was considered to be significantly enriched with differentially expressed genes.
Metabolic pathways enriched with differentially expressed genes. The KEGG database provides data and comments about metabolic pathways, allowing a visual representation of a gene's location in a metabolic pathway. The differentially expressed genes were compared to the KEGG database. According to the software calculation, get the p-value. If the p-value for a given metabolic pathway was less than 0.5, then the pathway was considered to be significantly enriched with differentially expressed genes.

Validation of the transcriptome database
After the total RNA was extracted from the skin tissue, the cDNA first chain was synthesized by reverse transcription using the FastQuant RT Kit (First-strand of cDNA synthesis kit, Tian-Gen Biotech(Beijing)CO.,LTD). The primer specificity was then detected.
Primer design and synthesis. Using the primer design program Primer Premier 5 and the sequence information in GenBank for AOC3, DUSP1, VCL, WNK1, IGFBP, and S100A7, fluorescence quantitative primers were obtained. The specificity of the primers was tested by Primer-BLAST. Primers for the internal gene sequences of goat GAPDH were obtained from a previous study by Chen Li[19]. All the primers were synthesized by Shanghai Sangon Biological Co., LTD. Table 1 shows the information for the primer pairs. PCR amplification products. Template cDNA and 2×Taq PCR MasterMix were used to synthesize seven genes up stream and downstream of the primers. The amplification system consisted of: 1 μL cDNA template, 1 μL forward primer (10 μM), 1 μL reverse primer (10 μM), 10 μL 2×Taq PCR MasterMix, and RNase-Free ddH2O added to make a total volume of 20 μL.
The response procedures were: 3 min at 94˚C; 30 cycles of 30 s at 94˚C, 30 s at 55˚C, and 1 min at 72˚C; and 5 min extension at 72˚C. The reaction product was stored at -20˚C.
Then, agarose gel electrophoresis and fluorescence quantitative PCR were performed. The 7500 Software v2.0.6 of the ABI7500 fluorescence quantitative PCR Software was used for data analysis. The 2 -ΔΔC method was used to calculate the relative expression of gene in triplicate. We used SPSS version 15.0 (SPSS, Inc., Chicago, IL, USA) for identify significant differences.
Screening of genes important for high-quality brush hair. The differentially expressed genes in the skin samples were determined from the high-quality hair and non high quality hair goats based on the standard of two folds change, and screened by the significant GO functional enrichment and metabolic pathway enrichment that might be important for traits related to high-quality brush hair.

RNA concentration and integrity testing
The result is shown in Fig 1. The RNA sample concentrations' result show as Table 2.

Transcriptome sequence assembly and gene functional annotation
Results of de novo assembly. The results of de novo assembly show as Table 3.
Gene functional annotation. The comparison to the NR, GO, KEGG, SwissProt, and egg-NOG databases resulted in hits to 32,812 genes in the NR database, 23,928 genes in the GO database, 12,703 genes in the KEGG database, 28,768 genes in the SwissProt database, and 29,337 genes in the eggNOG database.
Unigene eggNOG analysis. The eggNOG database had statistics for 29,337 of the genes in the skin samples divided into 25 functional categories. The most highly enriched category was Signal Transduction Mechanisms, followed by the Transcription and General Function categories (Fig 2).

Differential gene expression
Two hundred ninety-five genes were differentially expressed between the high hair quality and low hair quality goats. In the high hair quality goats, 132 genes were up-regulated, 62 genes were down-regulated, 101 genes were expressed specifically. Table 4 shows the top 20 up-regulated genes in the high hair quality goats.
GO functional enrichment of differentially expressed genes Table 5 shows the significant GO functional enrichment of the differentially expressed genes.
KEGG pathway enrichment with differentially expressed genes KEGG pathway analysis for significant enrichment with differentially expressed genes can predict involvement in biological function. The 295 differentially expressed genes were involved  Follow-up experiments were performed after the RNA quality was tested to meet the requirements for subsequent sequencing-library construction.
in 34 signaling pathways (Fig 3). Table 6 shows the KEGG pathways that were significantly enriched with differentially expressed genes. The KEGG pathway enrichment can show how the differentially expressed genes are involved in major functional and signal transduction pathways.

Important genes and pathways affecting type III hair growth
The top 20 differentially expressed genes that were up-regulated in the high-quality hair goats were conducted for the metabolic pathway analysis, of which three differentially expressed genes (MAP3K1, DUSP1 and DUSP6) were associated with the MAPK pathway, suggesting that the MAPK pathway plays an important role in the formation of type III hair. Furthermore, the MAPK signal pathway had totally six differentially expressed genes, with five up-regulated genes: MAP3K1, DUSP1, DUSP6, FOS and HSPA1, and one down-regulated gene: FLNA (Fig 4).

Fluorescence quantitative results
Agarose gel electrophoresis of PCR products. Agarose gel electrophoresis of the PCR products showed that all of the PCR products had a single band, indicating that the primer specificity, PCR reaction system, and PCR reaction program were satisfactory (Fig 5). Followup tests were therefore carried out.
Gene dissolution curve. After amplification, the dissolution curves of AOC3, DUSP1, GAPDH, IGFBP, VCL, WNK1, and S100A7 were analyzed (Fig 6). The dissolution curves each had a single peak, indicating that the specificity of the PCR products was good and that primer-dimers and non-specific PCR amplification did not occur. Relative gene expression. Assuming that the expression levels of VCL, WNK1, DUSP1, AOC3, IGFBP, and S100A7 were equal to 1 in the non high hair quality goats, the relative expression levels of those genes in the high hair quality goats were 5.67±0.25, 4.65±0.25, 2.92 ±0.18, 2.07±0.11, 0.78±0.16, and 0.66±0.13, respectively. The relative expression levels of VCL, Note: In the NCBI database, the warm-temperature-acclimation-related 65-kDa protein gene is not annotated in Capra hircus (goat). Using BLAST, we found the wap65-like gene (100049497) in Scombermorus niphonius (Japanese Spanish mackerel), which is 97% identical to the warm-temperatureacclimation-related 65-kDa protein gene(S1 File). Therefore, we used the gene and chromosome IDs of wap65-like in Table 4.
doi:10.1371/journal.pone.0169820.t004 Table 5. GO functional enrichment with differentially expressed genes. WNK1, DUSP1, and AOC3 were significantly higher in the high hair quality goats than in the non high hair quality goats (P<0.05). The relative expression levels of IGFBP and S100A7 were significantly lower in the high hair quality goats than in the non high hair quality goats (P<0.05). The relative expression trends were consistent with the results of transcription, indicating that the results of transcription were reliable, as shown in Fig 7.

Discussion
The RNA-Seq technique has been applied to many species and fields, such as Capra hircus [20], mice [17], and so on. It can quickly, comprehensively, and accurately detect the expression of all the genes in a sample over a certain period of time [21]. In this study, the transcriptomes of Yangtze River Delta White Goats with and without specific traits for high-quality brush hair were sequenced and analyzed. All of the Unigene annotation information for the respective samples was obtained. Analysis of the expression of all the genes in the samples identified 295 differentially expressed genes. In the goats with the high hair quality traits, 132 genes were up-regulated, 62 genes were down-regulated, and 101 genes were specifically  expressed. To verify the accuracy of the transcriptome database after high-throughput sequencing, most researchers have used real-time PCR technology [22][23][24][25]. Six differentially expressed genes, including four up-regulated genes and two down-regulated genes, were randomly selected in this study. The results of those six differentially expressed genes indicated that the sequencing results were reliable. The GO functional enrichment analysis of differentially expressed genes helps to describe the biological functions of genes of interest [26][27][28]. The GO analysis of the 295 differentially expressed genes resulted in 11 significantly enriched GO entries. In the Molecular Function The KEGG pathway analysis showed significant enrichment by the 295 differentially expressed genes in the Carbohydrate Metabolism, Energy Metabolism, and Terpenoid and Polyketone Compounds pathways, which accounted for 14% of the total differentially expressed genes. The differentially expressed genes accounted for 12% of the total differentially expressed doi:10.1371/journal.pone.0169820.g006 genes. Those correlations are higher when combined with GO analysis, the traits and metabolism, and the copy and translation of the genes in the Yangtze River Delta White Goats with the high-quality hair traits. It can be further speculated that the traits could be associated with cell proliferation and differentiation.
The morphology and structure of hair follicles, the parts of the skin that produce hair, are very complex. Hair follicle development is a finely-regulated process, and this process involves a large number of genes and pathways. Studies have shown that the MAPK signaling pathway is the main participants in hair follicle growth. The MAPK pathway plays an important role in biological functions such as cell proliferation and differentiation [29][30][31]. The MAPK signaling pathway also play an important role in the periodic development of hair follicles, which can promote the proliferation and differentiation of hair-follicle stem cells when it activated [32][33]. MAPK is one of the important signals mediating the cellular biological response system and is common in mammals. The MAPK pathway exerts its physiological and pathological regulation via the phosphorylation of substrate. The substrates of MAPK include transcription factors, intracellular protein kinases, and skeleton-related proteins, which highlights the diversity and complexity of MAPK functions [34]. At present, there are four known MAPK signaltransduction pathways in mammalian cells: EKR (extracellular-regulated protein kinase), JNK/SAPK (c-Jun N-terminal kinase/stress-activated protein kinase), p38, and ERK5/BMK1 [35].
The transcriptome sequencing results showed that three of the top 20 up-regulated differentially expressed genes, MAP3K1, DUSP6 and DUSP1, were involved in the MAPK signaling pathway, demonstrating that the MAPK signaling pathway likely plays an important role in the formation of type III hair in the Yangtze River Delta White Goat.
MAP3K1 was one of six differentially expressed genes in the MAPK signaling pathway. MEKK1 is the expressed protein of MAP3K1 and is an important node in the MAPK signaling pathway. That pathway can activate c-JUN, which enhances the transcriptional activity of transcription factor AP-1, thereby regulating cell migration and the formation of actin protein fibers [36]. The current research conducted differential proteomics analysis on high quality brush hair trait in Yangtze River Delta White Goat with the results showing that β-actin interacts with keratin to participate in the intracellular transport, cell shape, movement and division of hair follicle development related cells, which might play an important role in the process for promoting hair follicle development of type III hair [37]. Actin can affect the traits of type III hair, and MAP3K1 can regulate the formation of actin protein fibers; therefore, we believe that MAP3K1 is important in determining the traits of type III hair.
The dual-specificity phosphatase (DUSP) family proteins can regulate the activity of the MAPK pathway by dephosphorylation of key pathway proteins. This family includes DUSP1 and DUSP6. DUSP1 is also called MAPK phosphatase 1. It mainly catalyzes the hydrolysis of the phosphate group in the specific sequence TxY of the activated MAPK family members (p38, JNK, ERK) in the cell, thereby inhibiting the activity of the phosphoric acid group. The activation of MAPK can induce the expression of DUSP1. It is an important negative feedback regulation mechanism in the MAPK signaling pathway [38][39]. DUSP6 is also called MKP-3 (MAP kinase phosphatases 3). Some studies have shown that DUSP6 plays an important physiological regulatory function in the MAPK pathway [40][41]. DUSP6 can be combined with ERK2, and it can make ERK2 inactivation by the method of negative feedback control to regulate the activity of the MAPK pathway. DUSP1 and DUSP6 regulate the activity of the MAPK pathway, which is closely related to the development of hair follicles. For this reason, we suspect DUSP1 and DUSP6 may also be key genes that affect the quality of type III hair. At present, there have been many studies of DUSP1 and DUSP6 in cancer, but far fewer studies on the relationship between these proteins and the development of hair follicles. We will do more in-depth research on these genes in the future, and we hope these studies will enable a clearer understanding of their function in hair follicle development.
Gene function is not independent. Upregulation or downregulation of gene expression may affect the expression of several upstream or downstream genes, thus further altering gene expression patterns and the formation of some traits. We believe that specific expression patterns of these differentially expressed genes may affect the formation of type III hair traits.
To sum up, we think that MAP3K1, DUSP1, DUSP6 and MAPK pathway likely played an important role in the formation of the type III hair in the Yangtze River Delta White Goat.
Supporting Information S1 File. The result of blast in the warm-temperature-acclimation-related 65-kDa protein gene. (PDF)