Transcriptome Profile Analysis of Mechanisms of Black and White Plumage Determination in Black-Bone Chicken

In birds, feathers are skin appendages that play important roles in heat regulation, flight, and protection. Feathers have evolved to have different morphological and mechanical properties [1]. Plumage colour is a polygenetic trait involved in sexual selection, physical protection, geographical differentiation and speciation, and the evolution of polymorphisms [2-4]. Plumage colour is related to pigment distribution, content, and ratio [5]. Different plumage colours in birds are produced by different mechanisms. There are three distinct forms of plumage pigmentation, including melanin-based, carotenoid-based, and structurally based colours [6]. Melanin, which is endogenously synthesised in melanocytes and deposited in different tissues [7], is a ubiquitous component of plumage colouration in birds and animals. It comes in two main forms (eumelanin and phaeomelanin) in vertebrates, and yields black, brown, grey and earth-toned bird plumage colours [8]. Melanin pigments are synthesised by birds from the amino acids tyrosine, tryptophan, and phenylalanine, which are secreted by melanocytes and deposited in developing feathers [9]. Genetic abnormalities in the melanin pigment system can lead to a reduction in, or absence of, melanin synthesis, causing the animals to exhibit an albino phenotype (albinism) [10, 11]. A retroviral insertion in the TYR (tyrosinase) gene is associated with a recessive white mutation in chickens [10]. MC1R genes play a role in the rich diversity of plumage colour in chickens [12, 13]. Avian plumage colour is a complex trait involving major functional genes. However, the molecular mechanism of pigmentation is still largely unknown.

The Muchuan black-bone chicken is a local chicken breed in Muchuan County, Sichuan Province, China. It is characterised by an all-black body, including its feathers, cockscomb, skin, muscle, and claws. However, some individuals have black skin and white plumage; such that they differ from normal chickens only in their plumage colour. Black-bone chickens are hyperpigmented in various tissues, including the skin, muscle, ovary, and periosteum [14]. The special trait of hypermelanic pigmentation in the silky breed is due to dermal melanin and is independent of feather colour, as several varieties of this breed are described as having black skin and variable plumage colour: white, black, or gold [15]. In this study, a complete loss of melanin was observed in black-bone chicken feathers, whereas melanin pigment was synthesised normally in different chicken tissues, especially in the skin (Fig. 1). The variation in plumage colour (black and white) in black-bone chicken provides an excellent system for studying the molecular mechanism of melanin pigmentation during normal feather development. Study of the transcriptome profiles underlying the plumage colour of black-bone chicken, using high-throughput RNA deep-sequencing technology, may reveal a new mechanism of genetic control in the melanin plumage pigment system. This will provide a valuable avenue for research on melanin pigmentation in vertebrates.

A total of six female Muchuan black-boned chickens (52-week-old), including three black plumage and three white plumage chickens from a breeding farm in Muchuan County, were used in this study. All chickens were maintained at the same feeding and management levels. Feather bulbs (Fig. 1) from chickens with white (Group W: W1, W2 and W3) or black dorsal plumage (Group B: B1, B2 and B3) were collected and snap frozen in liquid nitrogen prior to use. All experiments involving animals were approved by the Leshan Normal University Animal Care and Use Committee. The methods were performed in accordance with the approved guidelines and regulations of the regional Animal Ethics Committee.

Total RNA was extracted using the Trizol reagent (Invitrogen, Life Technologies, Carlsbad, CA, USA) from the feather bulbs of chickens with black and white plumage according to the manufacturer’s instructions. The degradation and contamination of DNA was monitored on 1% agarose gels. The RNA purity and concentration were verified via the NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA) and Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies), respectively. RNA integrity was assessed via the RNA Nano 6000 Assay Kit in the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

A total amount of 3 µg RNA per sample was used to construct RNA sequencing libraries. Six sequencing libraries were constructed in total using the NEBNext Ultra RNA Library Prep Kit for the Illumina platform (NEB, Ipswich, MA, USA) according to the manufacturer’s recommendations, and index codes were added to attribute sequences for each sample. Briefly, we used poly-T oligo-attached magnetic beads to purify mRNA from total RNA. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext first strand synthesis reaction buffer (5×). The first cDNA and second strand cDNA were synthesised with random hexamer primer and M-MuLV reverse transcriptase (RNase H-), and DNA polymerase I and RNase H, respectively. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. NEBNext adaptors with hairpin loop structures were ligated to prepare for hybridisation after adenylation of the 3’ ends of DNA fragments. To preferentially select 150–200 bp-length cDNA fragments, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Prior to PCR, 3 µl of USER Enzyme (NEB) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min, followed by 5 min at 95°C. We then used Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer for PCR. Finally, PCR products were purified using the AMPure XP system and library quality was assessed using the Agilent Bioanalyzer 2100 system.

Following the manufacturer’s instructions, clustering of the index-coded samples was performed on a cBot cluster generation system using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina). Following cluster generation, library preparations were sequenced on an Illumina Hiseq X Ten platform and 150 bp paired-end reads were generated.

We processed raw data (raw reads) in fastq format via in-house Perl scripts. To obtain clean reads, reads containing adapter, reads containing ploy-N, and low quality reads were then removed from the raw data. We also calculated Q20 (1% error rate), Q30 (0.1% error rate), and GC content. All downstream analyses were performed based on high-quality clean data. The chicken reference genome (ftp://ftp.ensembl.org/pub/release-86/fasta/gallus_gallus/dna/Gallus_gallus.Gallus_gallus-5.0.dna.toplevel.fa.gz) and gene model annotation files (ftp://ftp.ensembl.org/pub/release-86/gtf/gallus_gallus/Gallus_gallus. Gallus_gallus-5.0.86.gtf.gz) were downloaded directly from the genome website. The reference genome index was built using Bowtie (ver. 2.2.3) and paired-end clean reads were aligned to the reference genome with TopHat (ver. 2.0.12) [16]. We selected TopHat as the mapping tool because it can generate a database of splice junctions based on the gene model annotation file, thus providing a better mapping result than other non-splice mapping tools.

HTSeq (ver. 0.6.1) was used to count the reads mapped to each gene. The fragments per kilobase of transcript per million mapped reads (FPKM) of each gene was calculated based on its length and the read count mapped to it. Differential expression analysis between black and white feather bulbs was performed using the DESeq R package (ver. 1.18.0). Using a model based on the negative binomial distribution, DESeq provided statistical routines to determine differential expression in the digital gene expression data. The resulting P values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P value < 0.05 found by DESeq were designated as differentially expressed [17, 18].

Gene Ontology (GO) enrichment analysis of DEGs was implemented using the GOseq R package, which corrects for gene length bias. GO terms with corrected P values < 0.05 were considered significantly enriched by DEGs.

KEGG is a database resource for understanding high-level functions and utilities of biological systems, 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-throughput experimental technologies (http://www.genome.jp/kegg/). We used the KOBAS software to test the statistical enrichment of DEGs in KEGG pathways [19].

The Cufflinks (ver. 2.1.1) Reference Annotation Based Transcript (RABT) assembly method was used to identify and construct both novel and known transcripts from TopHat alignment results. The alternative splicing (AS)events were classified to 12 basic types using the rMATS software (ver. 3.0.8). The number of AS events in each sample was estimated separately.

Picard-tools (ver. 1.96) and SAMtools (ver. 0.1.18) were used to sort and mark duplicated reads, and reorder the BAM alignment results of each sample. GATK2 (ver. 3.2) software was used to perform SNP and indel calling.

To validate DEG data, qRT-PCR was performed to determine the expression level of 12 randomly selected DEGs. Total RNA derived from the same samples used for RNA sequencing was reverse-transcribed to cDNA using the ProtoScript first-strand cDNA synthesis kit (NEB). The qRT-PCR primers were designed by Primer Premier 5 software (PREMIER Biosoft, Palo Alto, CA, USA) and are shown (for all online suppl. material, see www.karger.com/doi/10.1159/000489644) in Suppl. Table S1. qRT-PCR was conducted using SYBR Green Master Mix (Vazyme, Jiangsu, China) in the StepOne Plus Real-Time PCR system (Applied Biosystems, Foster City, CA, USA). The 20-μL PCR reaction consisted of 1 μL cDNA, 0.4 μL of each primer (10 μmol), 0.4 μL of ROX Reference Dye, 10 μL of SYBR Green Master Mix, and 7.8 μL of nuclease-free water. The amplification conditions consisted of an initial denaturation at 95°C for 5 min and 40 cycles of amplification (95°C for 10 s and 60°C for 30 s). A melt curve analysis was performed from 60°C to 95°C by reading the plate every 0.3°C. Detection was performed three times for each sample. Gene expression levels were calculated by the 2^–ΔΔCT method using β-actin as an internal control [20]. An unpaired Student’s t-test was used to evaluate significant differences between the two groups. All data are shown as means ± standard error (SE). P values < 0.05 were defined as significant and P values < 0.01 were highly significant.

TF analysis of DEGs was based on known and predicted TFs from the animal TFDB 2.0 database. For listed in species, we directly extracted the TF list from the database; otherwise, BLASTX was used to align the target gene sequences to the protein sequences of known TFs

A total of six cDNA library preparations (three from white dorsal plumage: W1, W2, and W3; and three from black dorsal plumage: B1, B2, and B3) were sequenced on an Illumina Hiseq X Ten platform {San Diego, CA, USA}. As shown in Table 1, the RNA sequencing (RNA-Seq) data sets from all libraries were highly consistent. The guanine–cytosine (GC) content of each library was more than 50%. Approximately 64% of the total reads mapped uniquely to the chicken reference genome (ftp://ftp.ensembl.org/pub/release-86/fasta/gallus_gallus/dna/Gallus_gallus.Gallus_gallus-5.0.dna.toplevel.fa.gz). The ratio of Q20 and Q30 for sequencing was greater than 95% and 89%, respectively. These differences in transcriptome profile were reflected in the correlations among samples (see online suppl. material, Suppl. Fig. S1). Pearson correlation coefficient values ranged from 0.94 to 0.97. These results indicate that the sequencing data for each library were of high quality, ensuring reliable results in further bioinformatics analysis.

The FPKM were computed using HTSeq software (http://www-huber.embl.de/users/anders/HTSeq/doc/overviewhtml) to quantify gene expression levels in RNA-Seq data analysis. According to this analysis, an average of 14, 690 expressed genes were detected in each sample (range: 14, 373–15, 146 genes). Gene expression level was divided into five categories (see online suppl. material, Suppl. Table S2): low expression (≤ 1 FPKM), low-to-moderate expression (1–3 FPKM), moderate expression (3–15 FPKM), moderate-to-high expression (15–60 FPKM) and high expression (≥ 60 FPKM). Genes with FPKM > 1 accounted for approximately 53.05% of the total annotated genes. Genes with high expression values (≥ 60 FPKM) accounted for approximately 5.51% of the total annotated genes. However, approximately 46.95 % of the expressed genes had low expression values (0 < FPKM ≤ 1). A total of 14, 531 and 15, 086 expression genes were detected in white (group W) and black (group B) feather bulbs, respectively, when the gene expression value of FPKM was larger than or equal to 1. In total, 14, 269 overlapped genes were filtered out from the total number of expression genes detected (see online suppl. material, Suppl. Fig. S2).

High-throughput RNA sequencing can reveal novel genes and their isoforms. In the present study, we compared the assembled transcripts with annotated genomic transcripts from reference sequences. A total of 1, 592 novel genes and 3, 549 novel isoforms were predicted from six chicken sequencing libraries: 806 novel genes were encoded on the sense strand, and 786 novel genes were encoded on the antisense strand.

To improve the accuracy of the gene annotation information in the current database, the boundaries of known genes were refined by alignment of known transcripts with assembled transcripts from the RNA-Seq data. In total, 11, 027 genes were structurally refined, including 5, 507 genes refined on the sense strand, and 5, 520 genes refined on the antisense strand. Detailed gene annotation information for the refined genes is provided in (see online suppl. material) Suppl. Table S3.

Alternative splicing (AS) is a key biological process to regulate gene expression. In the present study, the “junction counts” and “reads on target” methods within the rMATS software were used to analyse the AS events in each library. Five types of AS event were detected in the chicken feather bulbs, including skipped exon (SE), mutually exclusive exon (MXE), alternative 5' splice site (A5SS), alternative 3' splice site (A3SS) and retained intron (RI). A total of 8, 417 SE, 654 MXE, 102 A5SS, 188 A3SS, and 434 RI events were detected in six libraries using only the “junction counts” method, whereas 8, 437 SE, 655 MXE, 104 A5SS, 190 A3SS, and 447 RI events were identified using the “junction counts” and “reads on target” methods. The significantly different AS events between group W and group B, including 135 SE, 45 MXE, 7 A5SS, 6 A3SS, and 11 RI, were detected using only the “junction counts” method, whereas 135 SE, 49 MXE, 7 A5SS, 6 A3SS, and 16 RI events were identified by both methods. The numbers of upregulated and downregulated AS events for both detection methods are shown in Fig. 2.

Putative single nucleotide polymorphisms (SNPs) and insertion/deletion (indels) were detected in each library according to the alignment of the sequencing reads on the existing reference sequence. A total of 238, 872, 255, 256, 262, 764, 285, 737, 296, 996, and 270, 053 putative SNPs, and 23, 208, 24, 366, 25, 265, 26, 950, 27, 563, and 25, 364 putative indels were detected in W1, W2, W3, B1, B2, and B3, respectively.

We used DESeq software to detect DEGs at a significance level of P < 0.05. In total, 61 DEGs were identified between the feather bulbs from black and white chicken feathers (see online suppl. material, Suppl. Table S4). The log₂ fold change values ranged from –5.1569 to 5.9971 in feather bulbs samples. Upregulated DEGs included 47 genes, and downregulated DEGs included 14 genes. A total of 20 significantly expressed novel genes were identified between groups B and W that did not match any known genes in the database. To illustrate the significant differences in gene expression in six samples within the identified 61 DEGs, a heatmap and volcano plot are shown in Fig. 3.

Transcription factors (TFs) are key regulators of gene expression. We identified three TFs within newly identified DEGs from the current database, including HOXA9 (ENSGALG00000028983), GTF2IRD1 (ENSGALG00000001263), and OSR1 (ENSGALG00000016473). A novel gene (Novel00926) was identified in the FOXI2 family using BLASTX analysis.

To further assess the functional roles of 61 DEGs in plumage melanin pigmentation, gene ontology (GO) enrichment analysis was performed using the GOseq software based on Wallenius’ non-central hyper-geometric distribution. As shown in Fig. 4, one GO biological process term (skeletal system development) was significantly enriched in nine DEGs (seven upregulated genes and two downregulated genes). The most abundant GO terms in both black and white feather bulbs consisted of single-multicellular organismal processes, multicellular organismal processes, developmental processes, and single-organismal developmental processes. Other enriched terms for molecular function and cellular components contained 1–2 DEGs; further details are shown in (see online suppl. material) Suppl. Table S5.

We then conducted a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the DEGs. The results of this analysis showed that the DEGs were significantly involved in two pathways (P < 0.05): the transforming growth factor-beta (TGF- β) signalling pathway and ascorbate and aldarate metabolism (see online suppl. material, Suppl. Table S6).

To experimentally validate the expression levels of DEGs identified from the RNA-Seq data, quantitative real-time PCR (qRT PCR) was performed to detect the mRNA levels of 12 DEGs in feather bulbs. These genes included seven upregulated (HOXB9, CHRDL1, NTN1, BKJ, CHAC1, LAMA2, and GPX3) and five downregulated genes (IRX5, BMP5, RGN, RYR2, and STMN2). As shown in Fig. 5, the qPCR results confirmed our RNA-seq analysis results, indicating the same expression tendency in 12 genes obtained via both methods.

The development of avian plumage colouration is extremely complex. Plumage colouration can be categorised as pigment-based or structural [7]. Melanin is the most common component of colouration in animals; it is synthesised in melanocytes and deposited as granules in various organs [14]. Different pigment patterns form by modulating the presence, arrangement, or differentiation of melanocytes [21]. Therefore, melanin-based pigmentation is associated with a series of functional genes. In the present study, RNA sequencing was performed to explore the chicken transcriptome in feather bulbs of six black-bone chickens: three with black plumage and three with white plumage. We generated more than 300 million sequence reads, corresponding to 43.33 Gb of raw sequence data. In total, 61 DEGs were identified and classified into GO and KEGG categories. By comparing the two types of feather bulb, we found numerous novel genes, novel isoforms, AS events, SNPs, and indels, which may play crucial roles in melanin-based plumage pigmentation in chickens. Our findings provide novel and valuable insights into the molecular mechanism of melanin synthesis and deposition in avian plumage.

The number of DEGs was found to differ significantly between black and white feather bulbs. We discovered four interesting candidate genes for melanin pigmentation, all belonging to the homeobox family, that deserve further attention: HOXB9, HOXC8, HOXA9, and HOXC 9. The homeobox genes encode a highly conserved family of TFs that play an important role in development, regulating numerous processes including apoptosis, receptor signalling, differentiation, motility, and angiogenesis [22]. The function of Hox genes during the development of skin appendages, including hair follicles and feathers, has been reported in previous studies [23-28]. HOXB9 was identified as a Wnt target gene [29]. Wnt signalling is essential for organogenesis of the skin and its appendages, such as hairs, feathers, and scales [30]. Moreover, Wnt/β-catenin signalling plays important roles in melanocyte development and differentiation [31-33]. Chicken HOXC8 is reported to be expressed in dorsal dermal and epidermal cells during the first stage of feather morphogenesis [34].

We identified two glutathione (GSH) metabolism-related genes, ChaC glutathione-specific gamma-glutamylcyclotransferase 1 (CHAC1) and glutathione Peroxidase 3 (GPX3), as promising candidate genes for chicken feather melanin pigmentation. GSH, a tripeptide thiol found in virtually all animal cells, serves as an agent regulating the process of melanogenesis [35-37]. GSH levels are closely associated with the deposition of melanin in the skin of humans and other mammals [36, 38, 39]. The lowest levels of reduced GSH were found to be related to eumelanin-type pigmentation, whereas the highest were found in skin with phaeomelanin-producing melanocytes [40]. CHAC1 encodes a member of the gamma-glutamylcyclotransferase family of proteins, which catalyses cleavage of GSH into 5-oxoproline and a Cys-Gly dipeptide. Kumar et al. showed the overexpression of CHAC1-caused GSH depletion in mammalian cells [41]. Crawford et al. revealed that overexpression of CHAC1 led to a large decrease in total GSH, which was alleviated in a CHAC1 catalytic mutant in HEK293 cells [42]. Therefore, CHAC1 expression is associated with GSH metabolism; CHAC1 plays a critically important role in the process of melanogenesis. GPX3 belongs to the GSH peroxidase family that encodes protein and can catalyse GSH to form glutathione disulphide (GSSG). GSH-dependent peroxidase has been known to be involved in the reduction of hydrogen peroxide during both eumelanin and phaeomelanin synthesis [39, 43]. Levels of GSH, glutathione peroxidase, and glutathione reductase have been found to be directly related to the pigmentation of melanoma cells, suggesting that redox processes mediated by GSH play an active role in melanogenesis regulation [44]. We found that GPX3 expression level plays an important role in feather melanogenesis in chickens.

RNA-seq analysis was used to thoroughly examine the feather bulb transcriptome in the Muchuan black-boned chicken, which will help to elucidate the molecular mechanisms determining feather melanogenesis. From the results of pathway analysis, we identified two significant pathways related to plumage melanogenesis: the TGF-β signalling pathway, and ascorbate and aldarate metabolism. The TGF-β superfamily of paracrine and autocrine signalling molecules (activin, bone morphogenic proteins, TGF-βs, and decapentaplegic) is responsible for the regulation of many cellular behaviours, including proliferation, differentiation, migration, adhesion, and apoptosis [45-47]. TGF-β is known to be capable of suppressing the growth of normal human melanocytes [48-50]. Kishi et al. found that TGF-β participates in the regulation of chick retinal pigment epithelial cell proliferation and melanin synthesis [51]. TGF-β treatment markedly decreases the cell number of human primary melanocytes, whereas melanoma cells are less responsive to TGF-β-induced growth inhibition than normal melanocytes [52]. Kim et al. demonstrated that TGF-β1 decreases melanin synthesis via delayed extracellular signal-regulated kinase activation using a spontaneously immortalised mouse melanocyte cell line [53]. In this study, two genes were found to be involved in the TGF- β signalling pathway. Higher levels of BMP5 expression were detected in white feather bulbs, whereas PITX2 was expressed at a high level in black feather bulbs. These genes from the TGF- β signalling pathway may play an important role in chicken melanin synthesis. Only one gene, regucalcin (RGN), was found in the ascorbate and aldarate metabolism pathway. RGN is a calcium-binding protein that plays a pivotal role in maintaining intracellular calcium homeostasis [54]. Moisá et al. revealed that the ascorbate and aldarate metabolism pathway was involved in regulating longissimus muscle growth in Angus cattle [55]. However, we found no evidence from previous studies to support a relationship between melanogenesis and the ascorbate and aldarate metabolism pathway. RGN might be a promising candidate gene for further exploration of the role of RGN expression in feather melanogenesis.

To compare DEGs in melanin pigmentation between skin and plumage, eight well known melanogenesis-related genes (MITF, TYR, KIT, OCA2, ASIP, MCIR, KITLG, and TYRP1 [56, 57]) from previous researches in chickens were checked in black and white feather bulbs using RNA-seq and quantitative real-time polymerase chain reaction (qRT-PCR). RNA-seq analysis showed that except for TYR and MITF (P < 0.05), no significant difference in expression was detected between black and white chicken feather bulbs (P> 0.05). Expression of TYRP1 was not detected in black or white feather bulbs. qRT-PCR was performed to confirm the results obtained by RNA-seq (see online suppl. material, Suppl. Fig. S3). Using transcriptome profiling, Zhang et al. demonstrated that TYR, KIT, ASIP, and OCA2 were significantly differentially expressed between the black and white skins of black-bone chicken, whereas there was no difference in the expression of KITLG, MITF, or MC1R between the two skins. Based on these results, although feathers are elaborate skin appendages, the mechanism for melanogenesis between skin and feather in black-bone chicken is significantly different. Li et al. showed that the expression of TYR and MITF was associated with the formation of black and white plumage in ducks [58]. Melanin pigmentation is an important component of vertebrate colouration that is often under strong genetic control [3]. Many functional genes are involved in the melanogenesis pathway. In this study, we applied RNA sequencing to identify several melanin-related genes in chicken feather bulbs. We found major candidate genes for the melanogenesis and pigmentation processes in the chicken feather.

Our study provides a comprehensive overview of the transcriptome of black and white chicken feather bulbs. A set of differentially expressed candidate genes and two potential metabolic pathways are involved in black-bone chicken plumage melanogenesis: four homeobox genes (HOXB9, HOXC8, HOXA9, and HOXC 9), two GSH metabolism-related genes (CHAC1 and GPX3), and the TGF-β signalling pathway. Two known genes, TYR and MITF, were also demonstrated to play a role in melanin formation. Overall, our results provide a theoretical foundation for elucidating the molecular mechanisms of black and white plumage determination.

This study was supported by Scientific Research Project of Leshan Normal University (Z16020), and A Project Supported by Scientific Research Fund of Sichuan Provincial Education Department (17CZ0016). This work was also funded by National Natural Science Foundation of China (81502803) and Natural Science Research Program of Jiangsu Province (15KJB330005).

The authors declare that no conflict of interests exists.

S. Yu, G. Wang contributed equally to this work.