Skin Transcriptome Profiles Associated with Skin Color in Chickens

Nutritional and medicinal benefits have been attributed to the consumption of tissues from the black-boned chickens in oriental countries. Lueyang black-boned chicken is one of the native chicken breeds. However, some birds may instead have white or lighter skin, which directly causes economic losses every year. Previous studies of pigmentation have focused on a number of genes that may play important roles in coat color regulation. Illumina2000 sequencing technology was used to catalog the global gene expression profiles in the skin of the Lueyang chicken with white versus black skin. A total of 18,608 unigenes were assembled from the reads obtained from the skin of the white and black chickens. A total of 649 known genes were differentially expressed in the black versus white chickens, with 314 genes that were up regulated and 335 genes that were down-regulated, and a total of 162 novel genes were differentially expressed in the black versus white chickens, consisting of 73 genes that were up-regulated (including 4 highly expressed genes that were expressed exclusively in the skin of the black chickens) and 89 genes that were down-regulated. There were also a total of 8 known coat-color genes expressed in previous studies (ASIP, TYR, KIT, TYRP1, OCA2, KITLG, MITF and MC1R). In this study, 4 of which showed greater expression in the black chickens, and several were up-regulated, such as KIT, ASIP, TYR and OCA2. To our surprise, KITLG, MITF and MC1R showed no significant difference in expression between the black- and white-skinned chickens, and the expression of TYRP1 was not detected in either skin color. The expression of ASIP, TYR, KIT, TYRP1, OCA2, KITLG, MITF and MC1R was validated by real-time quantitative polymerase chain reaction (qPCR), and the results of the qPCR were consistent with the RNA-seq. This study provides several candidate genes that may be associated with the development of black versus white skin. More importantly, the fact that the MC1R gene showed no significant difference in expression between the black and white chickens is of particular interest for future studies that aim to elucidate its functional role in the regulation of skin color.


Introduction
The skin color of chickens is an important economic trait. Normally, there are three skin colors found in chickens: white, yellow and black. Skin color is the most direct marker whether the bird is black-bone chicken or not. In oriental countries, nutritional and medicinal benefits have been attributed to the consumption of tissues from the black-boned chickens for thousands of years. Therefore, skin color is a key trait that contributes to significant economic value in terms of poultry production. The Lueyang black-boned chicken is one of the native chicken breeds of Lueyang County in the Shaanxi Province of China. This bird is typically composed of eight characteristic black parts: feathers, wing tips, beak, cockscomb, skin, bones, legs and claws. However, some birds may instead have white or lighter skin, which directly affects the selective breeding of the Lueyang chicken population and causes economic losses every year. The presence of pigmented skin among Lueyang chickens is significantly tied to their economic value and the speed of breeding.
Pigmentation is a complex trait that depends on genetics and other factors, including the environment and certain drugs [1][2][3]. In vertebrates, 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 [4]. Many genes have been found to play well-known roles in pigmentation, based on previous genome-wide association scans (GWAS), and the analysis of these genes has identified many single nucleotide polymorphisms (SNPs), e.g., ASIP, OCA2, TYR, MC1R, KITLG, TYRP1, SLC24A4, MITF and KIT [5][6][7][8][9]. 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. For example, recent studies revealed that MC1R and ASIP are major genes involved in determining coat color in sheep [10,11], and TYRP1 is a strong candidate gene for color variation in Soay sheep [12]. An expression analysis was performed on 10 genes related to melanocyte development in Silky Fowl and White Leghorn embryos, via qRT-PCR, and a regulatory network for melanocyte development was constructed based on the expression data [13]. Emaresi et al. (2013) considered that color variation was likely to stem from differences in the expression levels of genes belonging to the melanocortin in the tawny owl [14]. Previous studies only focused on the gene expression patterns in animals with different coat color. However, in our study, even the chick had the same coat color, but the skin color was different. The aim of this work was to study the transcriptome profiles in the skin of chickens with black versus white skin using high-throughput RNA deep-sequencing technology, to investigate the different expression profiles of the genes involved in skin pigmentation, then look for the main differences between black and white skin colors in Lueyang chickens. This will enable us to understand the molecular mechanisms involved in skin pigmentation and provide a valuable theoretical basis for the selection of the black trait during the selective breeding of the Lueyang chick.

Materials and Methods
Ethics Statement completed at the poultry farm at Northwest Agricultural and Forestry University, Shaanxi, China. Ten healthy, 16-week-old white and black female Lueyang chickens (5 birds per color) were selected for the sample collection. A piece of skin (8 mm in diameter) from the left back was collected and immediately placed in liquid nitrogen. Total RNA was extracted from the sample using Trizol reagent (Takara, Dalian, China) according to the manufacturer's instructions. The RNA integrity was evaluated using gel electrophoresis, and the RNA purity was checked via the OD 260/280 ratio and the RIN value. RNA samples with a RIN value greater than 8.0 and an OD 260/280 ratio greater than 1.9 were selected for deep sequencing. According the result of sequencing, some colored gene expressions were validated using Real time quantitative polymerase chain reaction (qPCR). β-actin was used as housekeeping gene.

Library Generation and Sequencing
Three RNA samples from either the black or white skin samples were pooled following mRNA isolation. The isolated mRNA was fragmented, followed by first-strand cDNA synthesis using random hexamer primers. The second-strand cDNA was synthesized using buffer (Invitrogen, 20 μL), dNTPs (0.25 mM / μL), RNaseH (0.05 U / μL) and DNA polymerase I (0.5 U / μL). The short cDNA fragments were purified using the QIAQuick PCR extraction kit (LianChuan Sciences, Hangzhou, China). The fragment ends were repaired, A-tailed and ligated to sequencing adaptors. Suitably sized (350 ± 50 bp) fragments were selected following an agarose gel electrophoresis and used as templates for the PCR amplification to generate an RNA-seq library. The sequencing of the library was performed using an Illumina HiSeq 2000 (LianChuan Sciences, Hangzhou, China). The raw reads were cleaned by removing the adaptors and low quality reads prior to assembly. The unigene assembly was carried out using the short reads assembly program. All assembled unigenes were compared with the proteins in the non-redundant (nr) protein database, using BLAST software with a significance threshold of E-value < 10 -5 . Functional categorization by gene ontology (GO) terms was carried out according to molecular function, biological process and cellular component ontologies with an E-value threshold of 10 -5 . The pathway assignments were performed by sequence searches against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and using the BLASTX algorithm with an E-value threshold of 10 -5 .

Differential gene expression profiling
The expression abundance of each assembled transcript was measured through Fragments per Kilobase of exon model per Million mapped reads (FPKM) values. All reads were mapped onto the non-redundant set of transcripts to quantify the abundance of assembled transcripts. Bowtie was used for read mapping and applied for FPKM based expression measurement. The expressions of each reads between sample pairs (BS vs WS) were calculated using the numbers of reads with a specific match. Between the two samples, a minimum of a two-fold difference in log 2 expression were considered as expression differences.

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 MC1R, TYR, KIT, ASIP, TYRP1, OCA2, KITLG, MITF used for the qPCR can be found in Table 1. β-actin was used as housekeeping gene. Quantitative real-time PCR was performed in triplicate on the Stratagene iQ5 system. The 12.5 μL PCR reaction included 6.25 μL SYBR Premix Ex Taq II (TaKaRa, Dalian,China), 0.25 μL (10 p moL / μL) specific forward primer, 0.25 μL (10 p moL / μL) reverse primer, 0.5 μL ROX reference dye, 0.25 μL (10 ng / μL) diluted cDNA and 5.25 μL RNase free water. Cycling parameters were 95°C for 10 min, followed by 37 cycles of 95°C for 15 sec, 57°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.

Statistical Analysis
In this study, the skin from three females was used to prepare one pooled RNA sample for each group of white or black skin. Two cDNA libraries were then constructed to perform the Illu-mina2000 deep sequencing. The reference genome used in this study can be found at Gallusgallus.Galgal14.72.dna.toplevel.fa.gz. The significance level was |log2 (Fold change)| >1 with a p value < = 0.05. It ensures an accurate selection of the differentially expressed genes. The relative amount of mRNA expression of each gene (expressed as adjusted Ct value) was analyzed by the Stratagene iQ5 system. Adjusted cycle threshold (C (t)) values were calculated by following equation: CT value = 2 -ΔΔt . The analysis of variance was performed with SPASS Version 19.0 software. Differences were considered significant at P value < 0.05.

Illumina Draft Reads
In this study, we set up two libraries. One was of black skin, another was of white skin. The schematic of the Illumina2000 deep sequencing and analysis from the two libraries are shown in Table 2.

Differentially Expressed Genes
A comparison of the gene expression showed that a total of 649 unigenes were differentially expressed between WS and BS (|log2 (Fold change)| > 1, p value </ = 0.05), with 314 up-regulated genes and 335 down-regulated genes. A total of 162 significantly expressed novel genes were identified between WS and BS (WS vs. BS), consisting of 73 genes that were up-regulated (including 4 highly expressed genes that are expressed exclusively in the black-skinned chickens) and 89 genes that were down-regulated. Moreover, 506 unigenes were significantly differentially expressed in both WS and BS. In addition, 62 and 81 unigenes that were uniquely expressed in WS and BS, respectively. All of the differentially expressed genes are illustrated in Fig 1. M vs. A plots can be useful to determine any systematic bias that may be present between conditions (A = log2 (Fold change); M = -log10 (P value)). Therefore, we used the log2 (Fold change) value combined with the-log10 (p value) to assess every gene expressed between WS and BS. Volcano plots exploring the relationship between the fold change and the significance are shown in Fig 2. A heatmap was produced to explore the differences in gene expression. In a heatmap, different colors represent different expression levels, and a darker color represents greater expression. There are many genes with well-known roles in pigmentation that have been identified in previous genome-wide association scans (e.g., ASIP, OCA2, TYR, MC1R, KITLG, TYRP1, SLC24A4 and KIT). Several of these genes showed significant differences in expression between the white-and black-skin samples used in this study; TYR, OCA2, ASIP and KIT are shown in Fig 3. MC1R, KITLG, GRAP2, and IRF6 also showed no differences in expression. From the SLC family, SLC16A9 and SLC16A5 were significantly up-regulated (WS vs. BS). However,   contained triplicates of 10 test individuals (5 birds per color), one no-template-control, and a log10 dilution series. Samples were randomly assigned to PCR plates. The results indicated that there were no significant differences in the expression of MC1R, KITLG and MITF between the black and white skin samples (Fig 4), but the expression of TYR, KIT, OCA2 and ASIP was significant (P value < = 0.05); the results of the qPCR were consistent with the RNA-seq.

KEGG Analysis of the Pathways
A KEGG enrichment analysis showed that the differentially expressed genes are significantly involved in six predicted pathways (p value </ = 0.01), including steroid hormone biosynthesis, vascular smooth muscle contraction, phenylalanine metabolism, glycolysis / gluconeogenesis, calcium signaling pathway and melanogenesis. All of the pathways and the related information are described in Table 3. In this study, we paid more attention on the pathway of melanogenesis (Fig 5). A total of 40 differentially expressed genes are involved in the melanogenesis pathway, which plays an important role in animal pigmentation (Table 3). Three crucial genes were identified in this pathway: KIT, TYR, and ASIP, which were significantly up-regulated in BS (p  value </ = 0.01). In contrast, expressions of the MC1R and MITF gene were not significantly different between BS and WS (P value>0.05) (Fig 3).

GO enrichment Analysis
We analyzed the GO (gene ontology) enrichment of the identified genes. A total of 649 significantly different genes were mapped to the Gene Ontology database (http://www.geneontology. org/), and a total of 222 GO terms were obtained, including 77 significantly different terms (p value < 0.05). These terms were categorized into one of three categories: molecular function, biological process and cellular component (Fig 6). In the molecular function category, 10 significant unigenes are involved in calcium-ion binding (p value = 0.025). The canonical Wnt receptor signaling pathway is significantly different (p value = 0.0042) in the biological process category. In the cellular process category, the extracellular region is the most significant (p value = 0.0001), containing 10 significantly different unigenes.

Discussion
The mechanisms of melanogenesis are complex, and a number of recent publications have summarized the current understanding of the process of melanin production in the skin [15].
A number of studies, including GWAS, have identified numerous polymorphisms controlling human hair, eye and skin color, such as MC1R, OCA2, TYR, TYRP1, TYRP2, ASIP, IRF4, SLC24A4 and SLC45A2 [6][7][8]. A large number of single nucleotide polymorphisms (SNPs) have been identified, very few SNPs have been examined in relation to pigmentary phenotypes. It is well known that the rate-limiting enzyme for the formation of melanin is tyrosinase (TYR), which is responsible for several oxidative steps in the synthesis of melanin [16,17]. In our study, among the differentially expressed skin color genes, TYR showed the greatest level of differential expression in the black versus white chicken skin. This is consistent with studies of sheep coat color [18,19]. KIT is a receptor tyrosine kinase primarily composed of five extracellular immunoglobulin-like domains, a transmembrane domain and two intracellular tyrosine kinase domains. A new study showed that a novel splicing mutation in KIT results in piebaldism and auburn hair color in humans [20]. In regards to the important role of KIT in UVB-induced melanogenesis in the epidermis [21], more and more studies are focusing on this gene and the potential to lighten human skin color through the inhibition of KIT expression. In this study, the expression of the KIT gene was also significantly greater in the blackskinned compared to the white-skinned chickens. Black skin color is known to be due to an increase in the ratio of melanin relative to the white skin color. Previous studies revealed that MITF (microphthalmia-associated transcription factor) is involved in the development of melanocytes. The present study described the mutation at the MITF gene responsible for the plumage color both in the Japanese quail and gallus. The semi dominant B mutation results from a premature stop codon caused by a 2 bp deletion in exon11 of MITF [9]. The latest report indicated that there was alternative splicing of MITF gene in the skin of sheep [22]. It will provide new way to investigate the different color skin of chickens. Melanocortin-1 receptor (MC1R) is responsible for binding melanocyte stimulating hormone (MSH), which is expressed by stressed keratinocytes, and initiating the cascade of melanogenesis [23]. MC1R has classically been considered to play a role in and be the irreplaceable target involved in regulating mammalian skin pigmentation and hair color [24][25][26]. Despite the fact that more than 100 variants of the MC1R gene have been identified to date, the consequences of these variants on the physiological function of MC1R have been defined only partially [27]. Kinsler et al. (2012) hypothesized that if MC1R variants are associated with the growth of moles, the variants might also be associated with fetal growth in general [28]. It has been suggested that MC1R is likely to play an important role in the regulation of the immune system. Recent results have also shown that MC1R decreases inflammation both in vitro and in vivo and may be part of a therapeutic signaling pathway against inflammatory diseases [29]. All these findings expand on the known variations associated with the MC1R gene. More and more studies described that the MC1R gene had multiple pleiotropic effects that extend well beyond pigmentation. It can be summarized that MC1R does not stop at coloration. Instead, it has been proposed that this gene may also affect behavior, immune function, the nervous system and stress responses [30]. Therefore, many researchers have begun to pay more attention to the role of the MC1R gene in pigmentation, and recent studies have demonstrated new findings in regards to MC1R. Interestingly, the gene expression of MC1R, which has been identified in several studies as a gene related to pigmentation, showed no difference between the white and black skin samples. Three cosmetically important skin lightening agents, hydroquinone (HQ), kojic acid (KA), and niacinamide (NA) make up the bulk of the lightening ingredients in cosmetics, and an examination of these skin lighteners found that all of them led to marked increases in the expression of the TYR gene in the melanocytes but not MC1R [31]. Another recent study found that there is no association between MC1R polymorphisms and plumage coloration in the zebra finch [32], which is consistent with a study that was performed in leaf warblers. Another recently study also indicated that the coat color variation observed in rhesus macaques (variation from light to dark) is unlikely to be due to differences in the expression levels of the MC1R gene [33]. In a recent microarray experiment, no significant differences in the expression of MC1R were observed between the white skin and black spots of sheep. Instead, the researchers found many new differentially expressed genes that are likely to be involved in the formation of the black spots [34]. Similarly and to our surprise, MC1R showed no significant differences in expression between the black-skinned and white-skinned Lueyang chickens used in this study. MC1R is constitutively active so that there is no difference in expression between the two morphs. The specific function and role of MC1R in pigmentation still requires further research.
In a more recent paper, TYRP1 was found to show higher expression in the skin of blackcoated sheep versus white-coated sheep [19]. However, in our study, the expression of TYRP1 was not detected in either skin color. The results of our RNA-sequence and the subsequent validation were consistent, but these results still require additional 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 whiteand black-coated sheep [18], and a mutation in ASIP causes the black-and-tan pigmentation phenotype observed pigs [35]. Consistent with these previous studies, ASIP was also found to show significantly different expression levels in the black-skinned chickens versus the whiteskinned chickens used in this study. However, in the pathway of melanogenesis (Fig 5), It indicted that ASIP (ASP) gene was suppressed MC1R gene. In epidermal tissues, binding of ASIP to MC1R leads to the down-regulation of genes such as MITF and TYR. Finally, the production of pheomelanin would be reduced. In this study, the expression of the ASIP gene was also significantly greater in the black-skinned compared to the white-skinned chickens. So, it can suppress the expression of MC1R gene in black-skinned chickens. Whether it is the main reason for explaining MC1R gene had no significant difference between black-and white-skinned chickens or not, which need more investigations. Oculocutaneous albinism (OCA) is the most commonly inherited skin pigmentation disorder. Oculocutaneous albinism type 2 (OCA2), the most common OCA type, is caused by mutations in the OCA gene. The recent discovery of a novel OCA2 gene variant for melanoma risk in a familial melanoma pedigree provides evidence for the detrimental effects of OCA2 gene mutations on pigmentation, which supports the existing GWAS data regarding the relevance of the OCA2 gene in melanoma predisposition [36]. There are exceptions for both OCA2 and MC1R, which also harbor non-synonymous mutations, rs1800414 (His615Arg) in OCA2 and rs885479 (Arg163Gln) in MC1R, which have been associated with skin pigmentation in East Asian populations [37][38][39]. However, in sheep skin, the genes associated with OCA were expressed but the majority did not show any differential expression associated with the coat color [39]. OCA2 showed significantly different expression between the white-and black-skinned chickens used in this study, and the expression of OCA2 was up-regulated in black chick, which proving that OCA2 is closely related to skin color.
The GO and KEGG pathway analysis of the differentially expressed genes revealed that the majority were associated with cell functions and cell components. Of particular interest were the pathways related to pigmentation and melanogenesis. In this study, in terms of the molecular function category, 4 significant genes were found to be involved in melanogenesis, and 10 significant unigenes are involved in calcium-ion binding (p value = 0.025). In the biological process category, the canonical Wnt receptor signaling pathway showed a significant difference as well (p value = 0.0042). Both of these categories are related to pigmentation during the process of melanogenesis. In melanogenesis pathway, Wnt signaling pathway and Ca 2+ had important roles during melamin synthesis (Fig 5). All these results provide strong evidence that there is a significant difference in the level of melanogenesis between black-and white-skinned Lueyang chickens. However, further investigation is still needed to confirm the regulatory relationships of these genes.

Conclusion
In summary, this is an original report of the transcriptome analysis of the skin from chickens with white or black skin color. The present study described and revealed a set differentially expressed known and novel genes found in chicken skin that are potentially related to skin color and other physiological functions. The MC1R gene showed no difference in expression between chickens with black versus white skin color, and it is of particular interest for future studies to elucidate its functional roles in the regulation of skin color. The results provide a valuable theoretical basis for future breeding schemes for Lueyang chickens and a foundation for research regarding the manipulation of skin color.