Skin transcriptome profiles associated with coat color in goat

TL;DR

This study identifies genes associated with coat color in goats by analyzing skin transcriptomes from three different coat colors.

Contribution

The study provides new candidate genes and insights into the molecular mechanisms of coat color regulation in goats.

Findings

91 genes were differentially expressed between gray and white goat skin.

ASIP was only detected in white skin, while MC1R showed no significant expression differences.

qPCR validation confirmed RNA-seq results for most genes, except TYRP1.

Abstract

Capra hircus, an economically important livestock, plays an indispensable role in the world animal fiber industry. To identify additional genes that may play important roles in coat color regulation, Illumina/Solexa high throughput sequencing technology was used to catalog the global gene expression profiles in the skin of three different coat colors goat (Lubei white goat (white), Jining gray goat (gray) and Jianyang big ear goat (brown)). The RNA-Seq analysis generated 83174342, 70222592 and 52091212 clean reads in white skin, gray skin and brown skin, respectively, which provided abundant data for further analysis. A total of 91 genes were differentially expressed between the gray skin and white skin libraries, with 74 upregulated and 17 genes downregulated. Between the brown skin and white skin libraries, there were 23 upregulated genes and 44 downregulated genes, while there were 33 upregulated genes and 121 downregulated genes between the brown skin and gray skin libraries. To our surprise, MC1R, MITF, TYR, KIT and KITLG showed no significant difference in the skin of three different coat colors and the expression of ASIP was only detected in white skin and not in gray and brown skins. The expression of PMEL, TRPM1, DCT, TYRP1 and ELOVL3 was validated by real-time quantitative polymerase chain reaction (qPCR) and the results of the qPCR were consistent with the RNA-seq except the expression of TYRP1 between the gray skin and white skin libraries. This study provides several candidate genes that may be associated with the development of diferent coat colors goat skin. More importantly, the fact that the ASIP gene was only detected in the white skin and not in the other dark skins and the MC1R gene showed no significant difference in expression between the three different coat colors goat is of particular interest for future studies that aim to elucidate theirs functional role in the regulation of skin color. These results will expand our understanding of the complex molecular mechanisms of skin physiology and melanogenesis in goat and provide a foundation for future studies.

Figures

Classification of total raw reads pairs of different coat colors goat skin. After filtering the only adaptor sequences, containing N sequences and low quality sequences, the three RNA-Seq libraries still generated over 26 million clean reads pairs in each library, and the percentage of clean reads among raw tags in each library ranged from 94.07% to 94.71%.

Percent of reads mapped to genome regions. The reads mapped to genome regions including exon, intron and intergenic.

The numbers of DEGs between the three different coat colors goat skin. Between the gray skin and white skin libraries, there were 74 upregulated genes and 17 downregulated genes; Between the brown skin and white skin libraries, there were 23 upregulated genes and 44 downregulated genes, while there were 33 upregulated genes and 121 downregulated genes between the brown skin and gray skin libraries.

GO functional analysis of DEGs based on RNA-Seq data. The results were summarised the most enriched GO terms in three main categories: biological process, cellular component and molecular function.

The qPCR validation of differentially expressed genes in the three different coat colors goat skin. Abundance of target genes was normalized relative to abundance of GAPDH gene. Bars in each panel represent the mean ± standard error (n = 3).

Tables

Full text

The domestic goat (Capra aegagrus hircus) is a subspecies of goat domesticated from the wild goat of southwest Asia and Eastern Europe. The goat is a member of the family Bovidae and is closely related to the sheep as both are in the goat-antelope subfamily Caprinae. Goats have long been used for their milk, meat, hair and skins throughout the world (Wang et al. 2011a; Wang et al. 2011b; Xu et al. 2013).

Capra hircus, an economically important livestock, plays an indispensable role in the world animal fiber industry. Fiber diameter, length and color are key traits contributing to the economic value of goat and are determined by both genetics (Bunge et al. 1996; Lamoreux et al. 2001) and other factors, including the environment and certain drugs (Kidson and Fabian 1981; Dereure 2001; Jablonski and Chaplin 2010; Sturm and Duffy 2012). Factors that determine coat color in goat are becoming of increasing interest. White fleece holds greatest economic value due to its ability to be dyed to virtually any color, where as interest in natural colors is increasing due to the green revolution and consumer preference for natural products.

In adult animals, both hair and skin color depend on pigment produced by melanocytes at the base of the epithelium (Cieslak et al. 2011). Melanocytes in mammals and birds produce two types of melanin, black to brown eumelanin and yellow to reddish brown pheomelanin (Ito et al. 2000; Ito and Wakamatsu 2008). The basic coat colourations are defined by the ratio of the two pigments eumelanin and pheomelanin (Cieslak et al. 2011).

In animals, melanic coloration is often genetically determined and associated with various behavioral and physiological traits, suggesting that some genes involved in melanism may have pleiotropic effects (Ducrest et al. 2008). At present, a large number of genes have been found to play well-known roles in pigmentation and the analysis of these genes has identified many single nucleotide polymorphisms (SNPs), e.g., ASIP, MC1R, TYR, TYRP1, DCT, MITF, KIT, KITLG, OCA2, SLC24A4, PMEL (Brunberg et al. 2006; Gutierrez-Gil et al. 2007; Sulem et al. 2007; Deng et al. 2009; Nan et al. 2009; Duffy et al. 2010; Minvielle et al. 2010; Nicoloso et al. 2012; Hart et al. 2013; Fan et al. 2014; Li et al. 2014; Becker et al. 2015).

Coat color genes are good candidates for facilitation of trace ability of animal breeds. Several previous studies have paid significant attention to the coat color of animals and showed that this color is determined by the amount and type of melanin produced and released by the melanocytes present in the skin (Ito et al. 2000; Ito and Wakamatsu 2008). For example, recent researches have demonstrated that MC1R and ASIP are known to be major regulators of coat color in mice and MC1R and ASIP loci are functionally linked to undesirable coat color phenotypes in sheep (Vage et al. 1999; Slominski et al. 2004; Steingrimsson et al. 2006; Norris and Whan 2008). In addition, tyrosinase-related protein 1 (TYRP1) is a strong positional candidate gene for color variation in Soay sheep (Gratten et al. 2007). Recent studies have combined SNP analysis and gene expression profiling to dissect the basis for the piebald pigmentation phenotype in Merino sheep (Garcia-Gamez et al. 2011). It was reported that color variation was likely to stem from differences in the expression levels of genes belonging to the melanocortin in the tawny owl (Emaresi et al. 2013). Despite considerable knowledge of the genetic regulation of coat color in some animals, the molecular and cellular mechanisms regulating coat color in fiber-producing species, such as the goat, are not completely understood. This information is critical not only to enhanced basic understanding of regulation of melanogenesis, but also to the identification of novel pharmacological and molecular genetics approaches to regulate or select for coat color in fiber producing species.

RNA-Seq technology is a high-throughput sequencing platform allowing us to detect transcripts with low abundance, identify novel transcript units, and reveal their differential expression between different samples (Wang et al. 2009; Wilhelm and Landry 2009; Ozsolak and Milos 2011). To investigate the different expression profiles of the genes involved in goat skin pigmentation and gain insight into molecular mechanisms responsible for biochemistry of skin and fibers (including pigmentation) in animals producing hair, we investigated the transcriptome profiles in skin of goat of different coat colors using high throughput RNA deep sequencing. This will enable us to understand the molecular mechanisms involved in skin pigmentation and provide a valuable theoretical basis for the selection of the excellent natural colors trait.

MATERIALS AND METHODS

Ethics Statement

All of the animals were handled in strict accordance with good animal practices as defined by the relevant national and/or local animal welfare bodies. The experimental procedure was approved by the Animal Care and Use Committee of Hebei Normal University of Science and Technology, China and was performed in accordance with the animal welfare and ethics guidelines.

Goat skin sampling and total RNA extraction

Housing and care of goat and collection of skin samples for use in the described experiments were conducted in accordance with the International Guiding Principles for Biomedical Research Involving Animals (http://www.cioms.ch/frame 1985 texts of guidelines.htm). The animals were locally anaesthetized with hydrochloridum (1.5 ml of 3%, i.h.), following the approval provided by the Animal Hospital of Hebei Normal University of Science and Technology to decrease the animal suffering. 3 Lubei white goats, 3 Jining gray goats and 3 Jianyang big ear goats were selected for sample collection from the goat farm in Laiwu, Shandong province, China. A piece of skin (8 mm in diameter) from the leg was collected via punch skin biopsy under local anesthesia and immediately placed in liquid nitrogen. Total RNA from the sample was extracted using Trizol reagent (Invitrogen, USA) according to the manufacturer’s instructions. The RNA integrity was evaluated by gel electrophoresis and the RNA purity was checked by the ratio of OD260/OD280 and RIN value. RNA samples with RIN value greater than 7.5 and OD260/OD280 ratio greater than 1.75 were selected for goat sequencing.

Library generation and sequencing

Three RNA samples from white, gray and brown goat skin were pooled before mRNA isolation. Sequencing libraries were generated using NEBNext^®^ Ultra™ RNA Library Prep Kit for Illumina^®^ (NEB, USA) following manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer（5X）. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H^-^). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform and 100 bp paired-end reads were generated.

Reads mapping to the reference genome

Reference genome and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using Bowtie v2.0.6 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.9. We selected TopHat as the mapping tool for that TopHat can generate a database of splice junctions based on the gene model annotation file and thus a better mapping result than other non-splice mapping tools.

Differential expression analysis

The expression abundance of each assembled transcript was measured through Reads Per Kilobase of exon model per Million mapped reads (RPKM). Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the DEGSeq R package (1.12.0). The P values were adjusted using the Benjamini & Hochberg method. Corrected P-value of 0.005 and log2(Fold change) of 1 were set as the threshold for significantly differential expression.

GO and KEGG enrichment analysis of differentially expressed genes

Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with corrected Pvalue less than 0.05 were considered significantly enriched by differential expressed genes. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies (http://www.genome.jp/kegg/). We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways.

Real-time quantitative RT-PCR

The expression of these genes was quantified by qRT-PCR using QuantiTect SYBR Green RT-PCR (Qiagen, Waltham, MA). Information regarding the primers of PMEL, TRPM1, DCT, TYRP1, ELOVL3 used for the qPCR can be found in Table 1. GAPDH was used as housekeeping gene. Quantitative real-time PCR was performed in triplicate on the Stratagene iQ5 system. The 20 μL PCR reaction included 10 μL SYBR Premix Ex Taq II (TaKaRa, Dalian,China), 0.5 μL specific forward primer, 0.5 μL reverse primer, 0.4 μL ROX reference dye, 0.5 μL diluted cDNA and 8.6 μL RNase free water. Cycling parameters were 95°C for 4 min, followed by 40 cycles of 95°C for 15 sec, 56°C or 58°C for 30 sec and 72°C for 45 sec. Melting curve analyses were performed following amplifications. At the end of the cycles, melting temperatures of the PCR products was determined between 70°C and 90°C. The iQ5 software (Bio-Rad) was used for detection of fluorescent signals and melting temperature calculations. Quantification of selected mRNA transcript abundance was performed using the comparative threshold cycle (CT) method. The difference in abundance of mRNA for the genes was determined by analysis of variance.

Table 1.: Information regarding the specific primers used for the qPCR.

DATA AVAILABILITY

Table S1 contains top 30 highly expressed genes in goat skin. TableS2 covers differentially expressed genes between the gray and white skin. Table S3 covers differentially expressed genes between the brown and white skin. Table S4 covers differentially expressed genes between the brown and gray skin. Table S5 contains the list of GO categories for DEGs between different coat colors goat skin. Table S6 encompass KEGG pathway analysis for the gray and white skin differentially expressed genes. Table S7 encompass KEGG pathway analysis for the brown and white skin differentially expressed genes. Table S8 encompass KEGG pathway analysis for the brown and gray skin differentially expressed genes.

RESULTS

Massively Sequencing and Aligning to the Reference Genome

To maximise the coverage of the three different coat colors goat (Lubei white goat (white), Jining gray goat (gray) and Jianyang big ear goat (brown)) skin mRNA by RNA sequencing, RNA libraries were constructed by pooling RNA isolated from different coat colors individuals as a sample library. These RNA-Seq libraries generated over 54 million raw reads from each library. After filtering the only adaptor sequences, those containing N sequences and low quality sequences, the three RNA-Seq libraries still generated over 26 million paired-end clean reads in each library, and the percentage of paired-end clean reads among raw tags in each library ranged from 94.07% to 94.71% (Table 2 and Figure 1).

Table 2.: Quality Control of the Sequencing Data

Classification of total raw reads pairs of different coat colors goat skin. After filtering the only adaptor sequences, containing N sequences and low quality sequences, the three RNA-Seq libraries still generated over 26 million clean reads pairs in each library, and the percentage of clean reads among raw tags in each library ranged from 94.07% to 94.71%.

Of the total reads, more than 77% matched either to a unique or to multiple genomic locations; the remaining were unmatched (Table 3), because only reads aligning entirely inside exonic regions will be matched (reads from exon-exon junction regions will not be matched). Of the total mapped, more than 78% matched to exon and more than 11% matched to intron; the remaning were matched the intergenic (Figure 2).

Table 3.: Summary of read numbers based on the RNA-Seq data from different coat colors goat skin.

Percent of reads mapped to genome regions. The reads mapped to genome regions including exon, intron and intergenic.

Genes highly expressed in goat skin

The top 30 genes most highly expressed in goat skin included genes of the keratin family and ribosomal proteins (Table S1). The most highly expressed gene in all three different coat colors goat skin was keratin associated protein 3-1. The other highly expressed genes in the different coat color goat skins included, trichohyalin, transcript variant X1, keratin 71, keratin 5, keratin 25, eukaryotic translation elongation factor 1 alpha 1, keratin 14, transcript variant X1, keratin 27, keratin, type II microfibrillar, component 7C-like, collagen, type I, alpha 1, uncharacterized LOC102168701, ribosomal protein L35, secreted protein, acidic, cysteine-rich (osteonectin), ribosomal protein L13a.The list of most highly expressed genes also included 2 unknown genes.

Differentially expressed genes in different coat colors goat skin

To identify the differentially expressed genes (DEGs) of different coat colors goat skin, the differences in gene expression patterns were analysed for the pairs of the gray skin and white skin, the brown skin and white skin, the brown skin and gray skin. There were a total of 89 known genes and 2 novel genes identified as differentially expressed in gray skin versus white skin, of which 17 were down-regulated and 74 (including 2 novel genes) were up-regulated in skin from gray goat compared with skin from white goat (Figure 3 and Table S2). Between the brown and white libraries, 67 genes (including 1 novel gene) were differentially expressed, including 44 down-regulated (including 1 novel gene) and 23 up-regulated genes (Figure 3 and Table S3). When we compared the brown and gray libraries,154 DEGs (including 6 novel genes) were found, with 121 down-regulated (including 5 novel genes) and 33 up-regulated (including 1 novel gene) (Figure 3 and Table S4). This suggests that the differentiation of expressed genes between the brown skin and gray skin is larger than that between the gray skin and white skin, while a relatively smaller differentiation arises between the brown skin and white skin.

The numbers of DEGs between the three different coat colors goat skin. Between the gray skin and white skin libraries, there were 74 upregulated genes and 17 downregulated genes; Between the brown skin and white skin libraries, there were 23 upregulated genes and 44 downregulated genes, while there were 33 upregulated genes and 121 downregulated genes between the brown skin and gray skin libraries.

Functional Classification Analysis

For the GO analysis of the gray versus white skin, 1249, 349 and 607 differentially expressed genes were grouped in biological process, cellular component and molecular function, respectively. Most of the differentially expressed genes were classified into two GO categories (single-organism metabolic process, oxidoreductase activity and catalytic activity; Figure 4 and Table S5). The majority of the GO terms including pigmentation do not appear to be significantly enriched in the differentially expressed genes. For the GO analysis of brown versus white skin, 505, 326 and 391 differentially expressed genes were grouped in biological process, cellular component and molecular function, respectively. Most of the differentially expressed genes were classified into three GO categories (oxidation-reduction process, extracellular region and structural molecule activity; Figure 4 and Table S5). The majority of the GO terms including pigmentation do not appear to be significantly enriched in the differentially expressed genes. For the GO analysis of brown versus white skin, 2368, 761 and 1100 differentially expressed genes were grouped in Most of the differentially expressed genes were classified into two GO categories (oxidation-reduction process, single-organism metabolic process and oxidoreductase activity; Figure 4 and Table S5). The majority of the GO differentially expressed genes.

GO functional analysis of DEGs based on RNA-Seq data. The results were summarised the most enriched GO terms in three main categories: biological process, cellular component and molecular function.

In order to validate the transcriptome sequencing results, we selected 5 genes at random for real time PCR to determine their relative expression in different coat clor goat skin. These genes, identified as differentially expressed in dfferent coat color goat skin based on transcriptome sequencing analysis. The results of the qPCR were consistent with the RNA-seq except the expression of TYRP1 between the gray skin and white skin libraries (Figure 5).

The qPCR validation of differentially expressed genes in the three different coat colors goat skin. Abundance of target genes was normalized relative to abundance of GAPDH gene. Bars in each panel represent the mean ± standard error (n = 3).

KEGG pathway analysis

KEGG is a largely publicly available pathway-related database and is useful in searching for genes involved in metabolic or signal transduction pathways that were significantly enriched (Kanehisa et al. 2008). Of the 89 known genes differentially expressed in gray versus white goat skin, 40 had a specific KEGG pathway annotation. Of these KEGG pathway annotated genes, 3 were down-regulated in gray goat skin. These down-regulated genes are mainly involved in cysteine and methionine metabolism, melanogenesis, ribosome, metabolic pathways. Remaining KEGG pathway annotated genes were associated with 61 pathways including those functionally related to coat color in skin such as tyrosine metabolism (Table S6).

Between the brown and white goat skin, 23 DEGs with a KEGG pathway annotation were found and 14 were down-regulated in brown goat skin. These down-regulated genes were associated with 41 pathways including those functionally related to coat color in skin such as melanogenesis. Remaining KEGG pathway annotated genes were associated with 31 pathways including those functionally related to coat color in skin such as melanogenesis, tyrosine metabolism (Table S7).

Between the brown and gray goat skin, 72 DEGs with a KEGG pathway annotation were found and 63 were down-regulated in brown goat skin. These down-regulated genes were associated with 98 pathways including those functionally related to coat color in skin such as tyrosine metabolism. Remaining KEGG pathway annotated genes were associated with 29 pathways including those functionally related to coat color in skin such as melanogenesis, tyrosine metabolism (Table S8).

Differential expression of known coat color genes

Approximately 171 cloned genes involved in different pathways controlling coat color formation have been identified in the mouse (http://www.espcr.org/micemut/#cloned). Those known coat color genes are routinely classified into six general functions including: Melanocyte development, Components of melanosomes and their precursors, Melanosome construction/protein routing, Melanosome transport, Eumelanin and pheomelanin and Systemic effects (Bennett and Lamoreux 2003). Expression of a total of 128 of aforementioned coat color genes was detected in goat skin in present studies, and 2 genes showed higher expression in gray goat skin versus white goat skin and 1 gene showed higher expression in white goat skin, 4 genes showed higher expression in brown goat skin versus white goat skin and 1 gene showed higher expression in white goat skin, 4 genes showed higher expression in brown goat skin versus gray goat skin. Interestingly, the expression of ASIP (agouti signaling protein) were almost not detected in the gray skin and the brown skin, but there was high expression in the white goat skin for ASIP. PMEL, DCT, TRPM1 and TYRP1 were significantly up-regulated in brown goat skin (brown vs white, brown vs gray), but DCT, TRPM1 and TYRP1 showed no differences in expression between the gray and white skin samples and only PMEL was significantly up-regulated in gray goat skin versus white goat skin. ELOVL3 was significantly up-regulated in gray goat skin (gray vs white, brown vs gray). However, there were no significant differences in the expression of ELOVL3 between the brown and white skin samples. MC1R, MITF, TYR, KIT and KITLG showed no differences when comparing the three different coat color groups. The genes associated with oculocutaneous albinism (OCA) such as HPS1, HPS3, HPS4, HPS5 and HPS6 were expressed in goat skin but all of them did not show differential expression associated with coat color.

DISCUSSION

Mammalian coat color exhibits a wide range of shades and is dictated by melanin production in melanocytes (melanogenesis). Melanogenesis involves a complex molecular regulation (Slominski et al. 2004). In order to understand the molecular mechanisms of coat color formation, previous studies have reported the global gene expression profiles in skin of sheep and chicken with white versus black skin through Illumina sequencing technology (Fan et al. 2013; Zhang et al. 2015). Another study examined differences in gene expression associated with black spots infleece of Corriedale sheep using microarray technology (Penagaricano et al. 2012).

Our study offers new information related to gene expression profiles for three different coat colors goat skin. Our data analysis was based on the NCBI database of the goat reference genome (Capra hircus (assembly CHIR_1.0)).

To further investigate genes that may play important roles in goat skin, particularly in fiber/coat pigmentation, over 200 million transcriptome sequence reads were generated from different coat colors goat skin using high throughput RNA deep sequencing. From these reads there were 24,086 known genes identified as expressed in goat skin, among which 89 were differentially expressed in gray versus white goat skin, 66 were differentially expressed in brown versus white goat skin and 148 were differentially expressed in brown versus gray goat skin. It is acknowledged that study design was not optimal due to limited biological replication because single pooled samples (n=3 per coat color) were used in transcriptome sequencing analysis and the same three samples from different coat colors goat skin were used individually for quantitative real time PCR validation of the sequencing results. Despite such limitations, results have significantly enhanced understanding of goat skin transcriptome composition and potential differences in gene expression associated with coat color that are foundational to further study in the future.

Genes encoding for keratin family members and keratin associated proteins, ribosomal proteins were among the most highly expressed genes detected in goat skin. Of the top 30 highly expressed genes in goat skin, all 9 keratin family members and keratin associated proteins displayed no diffrenence between the three different coat colors goat skin. Hair keratins contain a much higher content of cysteine residues in their non-helical domains and thus form tougher and more durable structures via inter molecular disulfide bond formation (Mogosanu et al. 2014). Therefore, high expression of keratins is likely crucial for fleece strength. The ribosome is a central player in the translation system and its function is to decode the nucleotide sequence carried by the mRNA and convert it into anamino acid primary structure (Marshall et al. 2008). Abundant presence of ribosomal proteins in goat skin suggests the importance of high rates of protein translation in goat skin. In channel catfish skin, the expression of ribosomal proteins was high presumably due to higher levels of translational activities (Karsi et al. 2002; Patterson et al. 2003). It was also found that several collagen family members, particularly collagen, type I family members, were high expression in the three different coat colors goat skin. Collagen, in the form of elongated fibrils, is mostly found in fibrous tissues such as tendons, ligaments and skin (Zuber et al. 2015). Collagen has many medical uses in treating complications of the bones and skin. So high expression collagen may play important physiological function for skin. The expression of trichohyalin, transcript variant X1 was also high expression in the three different coat colors goat skin. Trichohyalin confers mechanical strength to the hair follicle inner root sheath and to other toughened epithelial tissues, such as the hard palate and filiform ridges of the tongue, by forming multiple complex crosslinks with itself and with other structural proteins (Steinert et al. 2003), therefore high expression trichohyalin plays important role in hair follicle development.

The GO and KEGG pathway analyses of differentially expressed genes revealed that most were associated with the function of single-organism metabolic process, oxidoreductase activity and oxidation-reduction process ontology categories. Of particular interest in our dataset were pathways related to pigmentation and melanogenesis. In this study, in terms of the molecular function category, 3 significant genes were found to be involved in melanogenesis and 3 significant genes were found to be involved in tyrosine metabolism. All these results provide strong evidence that there is a significant difference in the level of melanogenesis between the different coat colors goat skin. However, further investigation is still needed to confirm the regulatory relationships of these genes.

Melanocortin-1 receptor (MC1R) is responsible for binding melanocyte stimulating hormone (MSH), which is expressed by stressed keratinocytes, and initiating the cascade of melanogenesis (Solano et al. 2006). MC1R has classically been considered to play a role in and be the irreplaceable target involved in regulating mammalian skin pigmentation and hair color (Roberts et al. 2006; Schaffler et al. 2006; Lalueza-Fox et al. 2007). Recently, a study found that there is no association between MC1R polymorphisms and plumagecoloration in the zebra finch (Hoffman et al. 2014) which is consistent with another study that was performed in leaf warblers (Macdougall-Shackleton et al. 2003). Another recently study also indicated that the coat color variation observed in rhesusmacaques (variation from light to dark) is unlikely to be due to differences in the expression levels of the MC1R gene. In a recent microarray experiment, no significant differences in the expression of MC1R were observed between the white skin and black spots of sheep (Penagaricano et al. 2012). There were no significant differences in the expression of MC1R between the black-skinned and white-skinned Lueyang chickens (Zhang et al. 2015). Similarly and to our surprise, MC1R showed no significant differences in expression between the three different coat colors goat used in this study. The specific function and role of MC1R in pigmentation still requires further research.

The gene for agouti signaling protein (ASIP) is centrally involved in the expression of coat-color traits in animals. Previous studies have shown that the ASIP gene is responsible for the skin color of both white and black-coated sheep (Norris and Whan 2008) and the ASIP CNV were association with coat color in Tibetan sheep (Han et al. 2015). A mutation in ASIP causes the black-and-tan pigmentation phenotype observed pigs (Drogemuller et al. 2006), A missense mutation in the agouti signaling protein gene (ASIP) is associated with the no light points coat phenotype in donkeys (Abitbol et al. 2015). In our study, ASIP was found to show significantly different expression levels in the gray-skinned goats versus the white-skinned goats and the brown-skinned goats versus the white skinned goats, while ASIP was only detected in the white skin and almostly not detected in the gray skin and brown skin. The result that ASIP only was detected in the white skin and not in the dark skins seems to be inconsistent with the previous study, which has shown that the expression of the ASIP gene was also significantly greater in the black-skinned compared to the white-skinned chickens (Zhang et al. 2015). ASIP may play important role in the pigmentation of goat and further investigations are still needed to confirm the mechanism of ASIP.

TYRP1, one of themembers of the tyrosinase family, is a I type membrane bound protein that is expressed in both melanocytes and the retinal epithelium. TYRP1 is involved in the distal eumelanic pathway and plays a role in stabilizing TYR, which is the critical and rate-determining enzyme in melanogenesis (Kobayashi et al. 1998). There existed a significant association between coat color and TYRP1 in Soay sheep. In the free-living Soay sheep, coat color is either dark brown or light tawny color. The light phenotype is determined by homozygosity of a single recessive amino acid change (G→T transversion) at coding position 869 in the TYRP1 gene (Gratten et al. 2007). This is consistent with studies in domestic sheep, where light coat color is known to be due to a decrease in the ratio of eumelanin to pheomelanin, relative to black coat color (Aliev et al. 1990). In Bionda dell’Adamello goat, the genotype frequency distribution of a non-synonymous SNP (G1112A) suggested a possible role of TYRP1 in brown eumelanic goat coat colour (Nicoloso et al. 2012). The brown coat colour of Coppernecked goats is associated with a non-synonymous variant (p.Gly496Asp) at the TYRP1 locus (Becker et al. 2015). In our study, TYRP1 were significantly up-regulated in brown goat skin (brown vs white, brown vs gray) and TYRP1 showed no differences in expression between the gray and white skin samples using high throughput RNA deep sequencing, but the subsequent validation were inconsistent with our transcriptome sequencing data. The results of the qPCR showed that TYRP1 were significantly up-regulated in brown goat skin (brown vs white, brown vs gray) and had significantly higher expression in gary goat skin compared with white goat skin.These results still require additional research.

PMEL (SILV or PMEL17) is involved in melanosomal structure and acts as a scaffold in melanosomes by forming the proteolytic fibrillar matrix where melanin is deposited (Theos et al. 2005). It is a prominent member of the group of dilution genes. SILV mutations, which result in pigment dilution, have been discovered in many animals (Kerje et al. 2004; Brunberg et al. 2006; Clark et al. 2006; Gutierrez-Gil et al. 2007; Kuhn and Weikard 2007; Hellstrom et al. 2011). In our study, PMEL were significantly up-regulated in brown goat skin (brown vs white, brown vs gray) and had significantly higher expression in gary goat skin compared with white goat skin. The results of our RNA-seq and the subsequent validation were consistent.

DCT (TYRP2) is a melanogenic enzyme termed DOPAchrome tautomerase, which isomerizes the pigmented intermediate DOPAchrome to DHICA (5,6-dihydroxyindole-2-carboxylic acid) rather than to DHI (5,6-dihydroxyindole) (Tsukamoto et al. 1992). Mutations of DCT were associated with coat-color variation in mice (Jackson et al. 1992; Budd and Jackson 1995). The mice deficient in DCT showed a diluted coat color phenotype (Guyonneau et al. 2004). Mutations in dopachrome tautomerase (Dct) affect eumelanin/pheomelanin synthesis (Costin et al. 2005). Sheep with the GG genotype of DCT had significantly (P < 0.05) lower tyrosinase activity, alkali-soluble melanin content, and ratio of eumelanin: total melanin than sheep with GA and AA genotypes when measured across all investigated samples (Deng et al. 2009). DCT did not show expression in duck feather bulbs but expressed in retina (Li et al. 2012). Our RNA-seq and qPCR all showed that DCT were significantly up-regulated in brown goat skin (brown vs white, brown vs gray) and there were no significant differences in the expression of DCT between the gray and white skin samples.

TRPM1 is a member of the transient receptor potential (TRP) channel family. Although the potential effects of TRPM1 on pigmentation have not been fully elucidated, it has been reported to be a potent regulator of pigmentation (Oancea et al. 2009). Moreover, it is associated with a distinct leucistic white phenotype, the so-called leopard complex spotting, which is found in some horse breeds (Bellone et al. 2008; Wade et al. 2009). Interestingly, this mutation is linked with night blindness in homozygous Appaloosa horses. In homozygous leopard spotting horses, mRNA expression of TRPM1 is significantly downregulated in both skin and retina (Bellone et al. 2008; Bellone et al. 2010). This gene is also involved in night blindness in humans (Van Genderen et al. 2009; Koike et al. 2010). Our RNA-seq and qPCR all showed that TRPM1 were significantly up-regulated in brown goat skin (brown vs white, brown vs gray) and there were no significant differences in the expression of TRPM1 between the gray and white skin samples.

ELOVL3 encodes a protein that belongs to the GNS1/SUR4 family. Members of this family play a role in elongation of long chain fatty acids to provide precursors for synthesis of sphingolipids and ceramides. The Elovl3-ablated mice displayed a sparse hair coat (Westerberg et al. 2004). Our RNA-seq and qPCR all showed that ELOVL3 were significantly up-regulated in gray goat skin (gray vs white, brown vs gray) and there were no significant differences in the expression of ELOVL3 between the brown and white skin samples.

ASIP was only detected in the white skin and not in the dark skins, therefore ASIP may play important role in maintaining goat white coat color. PMEL,TYRP1, TRPM1 and DCT were significantly up-regulated in brown goat skin compare with the gray and white color skin, so they may play important role in the formation of goat brown coat color. In addition, ELOVL3 were significantly up-regulated in gray goat skin compare with the brown and white color skin and there were no significant differences in the expression of ELOVL3 between the brown and white skin samples. PMEL had significantly higher expression in gary goat skin compared with white goat skin and it was also found that TYRP1 had significantly higher expression in gary goat skin compared with white goat skin by qPCR. ELOVL3, PMEL and TYRP1 may be important for the goat gray coat color although the result of TYRP1 qPCR were inconsistent with RNA-seq. PMEL may be a key factor in the pigmentation of the three different coat colors. These results still require additional research.

CONCLUSIONS

In summary, this is an original report of the transcriptome analysis of the skin from three different coat colors goat. The present study described and revealed a set differentially expressed known and novel genes found in goat skin that are potentially related to skin color and other physiological functions. The MC1R gene showed no difference in expression between the three diferent coat colors and the ASIP gene was only detected in the white skin, which are of particular interest for future studies to elucidate its functional roles in the regulation of skin color. PMEL,TYRP1, TRPM1 and DCT may play important role in the formation of goat brown coat color and ELOVL3, PMEL and TYRP1 may be important for the goat gray coat color. Results are foundational for future studies to potentially manipulate coat color via pharmacological and genetic approaches.