Genome-Wide Association Study and Transcriptome Analysis Provide New Insights into the White/Red Earlobe Color Formation in Chicken

Chicken earlobe is a skin structure below the ear and covered without feather. As an external trait affected by both of natural selection and artificial selection, chicken earlobe color has accumulated lots of phenotypic variations (mainly including the red, the white, the predominantly white, the blue, the yellow and the black) in different breeds, being a typical breed characteristic in some breeds [1].

Compared with other earlobe color phenotypes in chicken, earlobe color with the white/ red has been relatively studied more, due to their predominance in population. Through crossing experiments in several breeds, Warren [2] found no inheritance mode could be concluded for white earlobe color, and it seemed to be sex-linked and incomplete dominance. Using 60K SNP chips, Wragg et al. [3] carried out a genome-wide association study (GWAS) and identified seven SNPs located on chromosomes 1, 2, 4, and Z that were associated with white/red earlobe color, further supporting the fact that the white earlobe was sex-linked and polygenic. Within a range of 2.38 Mb genomic regions on Chromosome Z, Nie et al. [4] also identified a total of 282 SNPs involving white/red earlobe color but without any other genomic region or locus on the left chromosomes, with 600K SNP chips used. Therefore, to the existing knowledge, the white/red earlobe color appears to be polygenic and sex-linked in some chicken breeds, and GWAS involving white/red earlobe color in more breeds will help further clarify its genetic mechanism. More importantly, until now we still know little about the molecular and histological mechanisms underlying the formation of white/red earlobe color.

Hematoxylin Eosin (HE) staining and melanin masson-fontana staining are two classical methods widely utilized in skin histology [5, 6]. Of the potential factors involving skin color in animals or human, the melanin, which is produced by melanocytes in epidermal basal layer, can play important roles in influencing skin pigmentation via hyperpigmentation or hypopigmentation [6, 7].

Qingyuan Partridge chicken (QYP) is a Chinese indigenous breed, of which the individuals with white earlobe color (Fig. 1A) and red earlobe color (Fig. 1B) coexist. In this study, we conducted histological sections on earlobe tissues of QYPs between 4 white earlobes and 4 red earlobes, aiming to investigate the potential histological difference between white and red earlobes. Meanwhile, to further clarify the inheritance mode and molecular genetic mechanism of white/red earlobe color, we performed a SLAF-GWAS (SLAF: Specific-Locus Amplified Fragment-sequencing) using 100 QYPs (30 males and 70 females) with white earlobe color and 100 QYPs (29 males and 71 females) with red earlobe color, and a transcriptome study of earlobe tissues between 3 114-day male QYPs with white earlobe color and 3 114-day male QYPs with red earlobe color.

The animal experiments in the present study were approved by the South China Agricultural University Institutional Animal Care and Use Committee (Guangzhou, People’s Republic of China) with the approval number SCAU#0011. All the experiments were handled in compliance with the guidelines of this committee.

A total of 200 300-day Qingyuan Partridge chickens (QYPs), which was obtained from Qingyuan FengXiangMaJi Co. LTD (Qingyuan, China), were used for GWAS to identify genomic regions associated with white/red earlobe color. QYP is a Chinese native breed which had four grades of earlobe color including red, predominately red, white and predominately white. Of the 200 chickens used here, 100 QYPs were with white earlobe (numbered from 1 to 100), whereas the other 100 QYPs were with red earlobe (numbered from 101 to 200). At 300 days of age, 0.3 mL blood samples were collected from each chicken for genomic DNA extraction.

According to our on-site observation, the white earlobe color phenotype just could be distinguished from other earlobe colors after 100 days in QYP. Therefore, we herein randomly selected 3 male 114-day QYPs exhibiting white earlobe color and 3 male 114-day QYPs exhibiting red earlobe for RNA-seq. In detail, the earlobe skins at a same position of the above 6 QYPs were collected and fast stored at -80 ℃ for RNA extraction.

Additionally, skin tissues from each 4 300-day male QYPs exhibiting red earlobe color (numbered with R_1, R_2, R_3 and R_4) and each 4 300-day male QYPs exhibiting white earlobe color (numbered with W_1, W_2, W_3 and W_4) were separately cut and stored in 4% paraformaldehyde for hematoxylin-eosin (HE) staining and melanin masson-fontana staining, to investigate the histological difference between red and white earlobe color groups.

HE staining and melanin masson-fontana staining was separately implemented as previously described by Lillie [8] and Kwon-Chung et al. [9].

Genomic DNA was isolated from whole blood using the NRBC Blood DNA Kit (Omega Bio-Tek, Norcross, GA, USA) following the manufacturer’s instruction. Total RNA of the earlobe skin was extracted with RNAiso Plus Kit (TaKaRa, Japan) used, in accordance with the manufacturer’s instruction. Both the DNA and RNA concentration were detected on a Nanodrop 2000 spectrophotometer at 260/280 nm ratio.

In this study, the red jungle fowl genome (Gallus_Gallus 4.0, ftp://ftp.ncbi.nlm.nih.gov/genomes/all/ GCA_000002315.2_Gallus_gallus-4.0) was chosen as the reference genome to perform enzyme digestion prediction with in-house related software. The fragments predictably digested by the enzyme that was selected for genome digestion should be in line with below rules, including a low proportion of that locating repetitive sequences, a uniform distribution across the genome to the greatest extent, a preferable identity between the fragment length and the experiment system [10] and the digested fragments number predicted meeting the expectation (i.e., 100, 000 SLAF). We herein selected restriction enzyme HaeIII to digest QYP DNA samples, predictably giving rise to 113, 202 SLAF, with the length ranging from 364-444 bp (for all online suppl. material, see www.karger.com/doi/10.1159/000489361, Table S1; Fig. S1). After digested by HaeIII, the genomic DNA fragments (SLAF) from each sample was added with A-tailing at 3’ end, further ligated with Dual-index adaptor, and then amplified by PCR. After that, SLAF from each sample was purified and then mixed for further target fragment selection via gel electrophoresis and extraction, thereby obtaining the SLAF library for high-sequencing. Ultimately, SLAF libraries from each sample were sequenced on Illumina HiSeq 2500 platform, according to Illumina instructions (operated by Biomarker Biological Technology Co., Ltd., Beijing, China). We herein developed a total of 622, 443 SLAF markers, including 424, 810 polymorphic ones (see online suppl. material, Table S2), and their corresponding average sequencing depth reached 10.03 × (see online suppl. material, Table S3).

Raw data was first processed using in-house Perl scripts, to remove reads containing adapter, reads containing over 10 % of ploy-N and low-quality reads (> 50 % of bases whose Phred scores were < 5 %), thereby obtaining clean reads. To enable the dependability of the reads for further analysis, the double-end 100 bp of the clean reads were chosen for subsequent analysis. Further, the filtered clean reads were assigned to map the reference genome Gallus Gallus 5.0 by the method of BWA [11]. With the use of GATK [12] and SAMtools [13], Li et al., 2009, the SNPs that were both identified by above tools were defined as the true SNPs, resulting in 2, 172, 483 population SNPs. After removing the integrity less than 80% and the minor allele frequency (MAF) less than 5% of each population SNP, a total of 230, 546 qualified SNPs, basically uniformly distributed on genome (see online suppl. material, Fig. S2), was kept for GWAS (see online suppl. material, Table S4).

The obtained RNA was treated with RNase free-DNaseI (Promega, USA) to remove any contaminated genomic DNA. For each sample, a total of 3 μg RNA was used for sequencing library preparation by NEBNext® Ultra^TM RNA Library Prep Kit for Illumina® (NEB, USA), according to the manufacturer's instruction. In order to select cDNA fragments of 150∼200 bp in length, the library was purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was added into size-selected, adaptor-ligated cDNA for PCR at 37 °C for 15 min followed by 5 min at 95 °C. At last, PCR products were purified and the library quality was assessed on the Agilent Bioanalyzer 2100 system. The sample clustering was carrried out on a cBot Cluster Generation System with the use of TruSeq PE Cluster Kit v3-cBot-HS (Illumina, USA) according to the manufacturer's instruction. After cluster generation, the library was sequenced on an Illumina Hiseq 2500 platform. Eventually, 150 bp paired-end reads were generated.

Those 230, 546 qualified population SNPs were first performed Principle Component Analysis (PCA) to determine whether stratification existed between white earlobe group and red earlobe group, using EIGENSOFT software [14]. For GWAS, association of each qualified SNP with earlobe color was performed by the software package PLINK (v1.09), using the case-control association analysis [15]. Chicken with white earlobe color were set as the case group, while birds with red earlobe color was the control one. Bonferroni method at α ≤ 0.05 and 1 level were separately used to determine the significant and suggestive threshold (i.e., 0.05/230, 546 and 1/230, 546) [16]. Genome-wide P values of association analysis was generated and exhibited via Manhattan plot using an in-house Perl script. Moreover, associated genes involving white/red earlobe color were identified within the upstream and downstream 50 Kb genomic regions of the significant SNPs.

For RNA-seq, paired-end clean reads were aligned to the chicken reference genome (ftp://ftp.ensembl.org/pub/release-88/fasta/gallus_gallus/) using TopHat v2.0.12 [17, 18] after the reference genome index built with Bowtie v2.2.3 [19, 20]. HTSeq v0.6.1 was used to count the read numbers mapped to each gene. Then, gene expression levels were calculated by the number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) [21]. Differential expression analysis of two groups (3 biological replicates for each) was performed using the DESeq R package (1.18.0) [22]. Genes that exhibited more than 1.5-fold changes between two groups and an adjusted P-value < 0.05 were assigned as the differentially expressed. Gene functional analyses, including GO, KEGG and HPO, were implemented by g:Profiler (http://biit.cs.ut.ee/gprofiler/index.cgi), with corrected P-value < 0.05 considered significant. Protein-Protein Interaction network analysis (PPI) was conducted by STRING (v10.5) (https://string-db.org/cgi/input. pl?input_page_show_search=on&UserId=pkHSI3fUSakh&sessionId=NjPXVsCoRQ0B) [23].

Additionally, Haploview 4.2 was utilized for Linkage Disequilibrium (LD) analysis, single SNP and haplotype association tests [24], with association significance tested by 10, 000 permutations. And the association test model was case (white earlobe individuals) - control (red earlobe individuals) model.

The RNA-seq and SLAF-GWAS data has been uploaded to Sequence Read Archive (SRA) database with the accession number GSE101622 and GSE110145, respectively.

To investigate the potential histological mechanism underlying the white/red earlobe color in QYP, we separately performed HE staining and melanin masson-fontana staining on white earlobe tissues from 4 male 350-day QYP and red earlobe tissues from 4 male 350-day QYP. Through HE staining, in the dermis, we could observe the major difference between the white group and the red group was the density of blood or capillary vessels in dermal papillary layer (Fig. 2). In the dermal papillary layer of red earlobes, much denser blood vessels occurred compared to that in white earlobes. Whereas in the epidermis, no significant difference in cell numbers or types could be observed in cuticle, granular layer, prickle cell layer and basal layer (stratum lucidum could be observed in neither of them). Furthermore, through melanin masson-fontana staining, there were no obvious melanins or melanocytes that could be observed from epidermal basal layer or prickle cell layer (see online suppl. material, Fig. S3).

PCA result showed that no obvious stratification existed between white earlobe and red earlobe groups (see online suppl. material, Fig. S4), thus those 230, 546 qualified SNPs from 200 QYPs were wholly aggregated for GWAS to investigate the genomic regions or genes involving white/red earlobe color. As a result, we identified a total of 44 SNPs that was significantly associated with the white/ red earlobe color based on the Manhattan plot (Fig. 3). These SNPs spread across 12 chromosomes, including chromosome 1, 2, 3, 4, 9, 10, 11, 13, 19, 20, 23 and Z (Table 1). And the most significantly associated SNPs, rs316568372 and rs316893415, were from chromosome 9, while the chromosome Z had the largest number of significantly associated SNPs. Focusing on the upstream and downstream 50 Kb genomic regions of above significantly associated SNPs, we further identified 73 genes (Table 1) potentially concerning white/red earlobe color. Of them, tumor protein 63 gene (TP63) was close to the two most significant SNPs, rs316568372 and rs316893415; phosphatidylinositol-4, 5-bisphosphate 3-kinase catalytic subunit beta gene (PIK3CB) was adjacent to the three significant SNPs of rs14667843, rs14667813 and rs315019040; beta-1, 4-galactosyltransferase 1 gene (B4GALT1) was close to the two significant SNPs of rs14780309 and rs313777363. GO analysis on these 73 significantly associated genes showed PIK3CB, B4GALT1 and TP63 were closely related to anatomical structure formation involved in morphogenesis, and particularly PIK3CB and B4GALT1 were involved in angiogenesis, which could directly account for the distinct histological observations between white earlobes and red earlobes.

From RNA-Seq, we obtained 46, 208, 932, 44, 992, 438, 43, 529, 604, 50, 294, 718, 48, 577, 764, and 47, 756, 246 raw reads for white_1, white_2, white_3, red_1, red_2, and red_3, respectively, giving rise to more than 6.26 Gb clean bases for each sample (see online suppl. material, Table S5). Gene expression comparison between red group and white group (Red_ vs_White) showed that there were 81 differentially expressed genes (DEGs) (fold changes ≥ 1.5; adjusted P-value < 0.05) between the two samples (see online suppl. material, Table S6), in which 27 were up-regulated and 54 were down-regulated in Reds compared to Whites. Further functional analysis on DEGs showed that, no significant KEGG pathways or GO terms could be clustered, but a number of genes including TP63 etc, could be clustered into HPO terms of Thickened skin, Epidermal thickening, Hyperkeratosis, Dry skin and Alopecia (Table 2).

To elucidate the validity of the DEGs identified by RNA-seq, we herein selected 7 DEGs, including TP63, AQP9, SLC38A4, KRT5, LECT1, NOX1 and VILL, to perform qPCR, using another 6 114-day male QYPs (including 3 white earlobes and 3 red earlobes), normalized with β-actin (Primers information was presented in see online suppl. material, Table S7). QPCR results showed that the expression patterns of the above genes were basically concordant with that of RNA-seq except LECT1, suggesting a high reliability of RNA-seq results (Fig. 4).

We compared the DEGs obtained from RNA-seq here with the significantly associated genes of white/red earlobe color identified in our GWAS and previous concerned GWAS conducted by Wragg et al. [3] and Nie et al. [4], respectively (genes identified by previous concerned GWAS not shown here), finding that just the TP63 identified by our GWAS significantly differentially expressed. Further, we integrated the DEGs here and the significantly associated genes identified by the present GWAS, to perform PPI analysis, and thereby investigated the pivotal genes involving white/red earlobe color formation. Consequently, a Protein-Protein Interaction network comprising 12 node genes were identified (P-value: 0.0867), in which PIK3CB and Cadherin 1 (CDH1) were the hub node gene (see online suppl. material, Fig. S5), suggesting potentially important roles in white/red earlobe color formation played by them. In addition, no repetitive candidate white/red earlobe color-related genes were found among present GWAS results, Wragg et al. [3] and Nie et al. [4].

TP63 as the candidate gene identified by present GWAS, it was also differentially expressed in RNA-seq, highlighting its potential roles in chicken white/red earlobe color formation. Therefore, we analyzed the mutation (including 5` and 3` UTR) of TP63 transcript between the white earlobe and red earlobe individuals, thereby investigating whether there was a specific SNP that could monitor the differential expression pattern between them or associated with white/red earlobe color. But no such a mutation could be observed (data not shown here). Further, genotyping on proximal 5` flanking region of TP63 was performed with the 200 GWAS population used. Within the region from Chr.9_14, 189, 720 to Chr.9_ 14, 191, 557 (primers P2 and P3 used; (see online suppl. material) in Table S7), a total of 10 SNPs was identified, namely rs317138318 (P2-1), rs317582344 (P2-2), rs317395129 (P2-3), rs315786018 (P2-4), rs314583828 (P2-5), rs740487702 (P2-6), Chr.9_14190820 (P3-1), Chr.9_14190824 (P3-2), Chr.9_14190829 (P3-3) and Chr.9_14190876 (P3-4). LD analysis showed that there were 3 blocks across these 10 SNPs (see online suppl. material, Fig. S6), including Block 1, Block 2 and Block 3. Of them, Block 1 consists of 4 SNP sites (P2-1, P2-2, P2-3 and P2-4), including haplotypes C-G-A-T, T-A-G-A, T-A-A-T and C-A-A-T; Block 2 consists of 2 SNP sites (P2-5 and P2-6), including haplotypes A-G, G-G and A-A; and Block 3 consists of 2 SNP sites (P3-1 and P3-3), including haplotypes C-T and T-C. Single SNP and haplotype association tests on above 10 SNPs and 3 blocks indicated that, allele T of P2-4 (adjusted P value=0.04), allele A of P2-6 (adjusted P value=0.04), haplotype T-A-G-A in Block 1 (adjusted P value=0.04) and haplotype A-A in Block 2 (adjusted P value=0.04) were significantly correlated with white earlobe color (Table 3).

We herein utilized SLAF-GWAS to identify the potential candidate genes related to white/ red earlobe color in 200 QYPs. As a result, we identified a total of 44 genomic regions from Chromosome 1, 2, 3, 4, 9, 10, 11, 13, 19, 20, 23 and Z, corresponding to 73 genes within 50 Kb upstream or downstream genomic regions of those SNP markers. Our results did support the fact that the white/red earlobe color was polygenic and sex-linked, which was consistent with previous studies conducted by Wragg et al. [3] and Nie et al. [4]. Importantly, the significant SNP markers or genomic regions involving white/red earlobe color identified here were absolutely distinct from those identified by Wragg et al. [3] and Nie et al. [4], strongly indicating the genetic heterogeneity of the white/red earlobe color in chicken. A caveat was that the big difference between what we observed and previous concerned findings might be also attributable to the methods used for genotyping. In the studies of Wragg et al [3]. and Nie et al. [4], 60K and 600K SNP chips were used, respectively. But compared to SNP chips, the method of SLAF we used here was based on enzyme digestion, inevitably bringing a non-uniform SNP coverage on genome, thus further inducing some potential significant SNPs being omitted or unraveled. Utilizing the whole-genome re-sequencing for genotyping in white/red earlobe color, would help clarify the ambiguity.

Through a transcriptome study of earlobe tissues between 3 114-day male QYPs with white earlobe color and 3 114-day male QYPs with red earlobe color, we identified a total of 81 DEGs, and some of them (e.g., TP63) could be significantly enriched in the HPO terms like Thickened skin. In addition, no DEGs involving the metabolism of melanin and purine compounds were identified, weakening the hypothesis that hypopigmentation or deposition of compounds made up of purine base [25] contributed to white earlobe formation. Histological section analysis between the white earlobes and red earlobes showed that, there were many more blood vessels in the dermal papillary layer of red earlobes, and no obvious melanins could be detected in both the white ones and the red ones, also suggesting melanins was not the main cause for white earlobe. Importantly, we found the individuals with white earlobe color had thicker earlobe skins than the ones with red earlobe in QYP, through touching by hand. Therefore, thicker skins and fewer blood vessels in dermal papillary layer might have major effects on the white earlobe formation, as they could act together to pale the redness of skin.

PIK3CB, B4GALT1 and TP63 as the three genes identified at the present GWAS, they were the first time to be reported in chicken white/red earlobe formation. Particularly, the PIK3CB and B4GALT1, are both closely associated with angiogenesis [26-29], which can well account for the quantity distinction of blood vessels from the histological observations between white earlobes and red earlobes. Besides, TP63 as a homologue of the tumor-suppressor P53 [30], its function is closely associated with epithelial development and morphogenesis [31, 32]. TP63 is also an active gene involving tumorgenesis [33]. In squamous cell carcinoma (SCC), ΔNp63, the predominant isoform of TP63, its over-expression in vitro and in vivo could up-regulate human beta-defensin (HβD), thus promoting tumor angiogenesis in SCC microenvironment [34]. Importantly, RNA-seq study between the white earlobes and the red earlobes indicated that, TP63 was significantly down-regulated in white earlobes compared to that in red ones. Although we didn't identify any tag mutation that could monitor the expression pattern between white earlobes and red earlobes or associated with white/red earlobe color in TP63 transcript, we identified 2 alleles (i.e., allele T of P2-4 and allele A of P2-6) and 2 haplotypes (i.e., T-A-G-A in Block 1 and A-A in Block 2) were significantly correlated with white earlobe formation in the proximal 5` flanking region of TP63, indirectly supporting the potentially important roles in white earlobe formation played by TP63. Therefore, it was reasonable for us to propose that down-regulated TP63 in white earlobes possibly led to thicken the skin and decrease the blood vessel numbers, thereby contributing to the white earlobe formation via paling the redness of skin in QYP. In addition to TP63, CDH1 as the hub gene identified by PPI, it is an anti-angiogenesis gene [35], and its ablation with Pten in mouse uterus can accelerate angiogenesis [36]. Although, in this study, CDH1 was down-regulated in white earlobes compared with in red earlobes, which was inconsistent with expected, it could be a candidate gene closely associated with white/ red earlobe color formation to be further explored.

In conclusion, we found the major histological difference between the white earlobes and the red ones was the blood vessel numbers in dermal papillary layer. And we identified a total of 44 SNPs from Chromosome 1, 2, 3, 4, 9, 10, 11, 13, 19, 20, 23 and Z, that was significantly associated with the chicken white/red earlobe color from GWAS, corresponding to 73 associated genes. Of them, PIK3CB, B4GALT1 and TP63, as the genes significantly associated with the white/red earlobe color, are closely associated with angiogenesis. Our results do also support the fact that the white/red earlobe color is polygenic and sex-linked in QYP. Additionally, through contrast of RNA-seq profiles between the white group and the red group, we identified a total of 81 differentially expressed genes (DEGs), of which a number of DEGs (e.g., TP63) could be significantly clustered into the HPO terms like Thickened Skin. In particular, TP63 as the significantly associated gene, it was also found to be differentially expressed between the white earlobes and the red ones. Combined with the existing knowledge of TP63 in epithelial development and tumor angiogenesis, we deduce that down-regulated TP63 in white earlobes may lead to thicken the skin and decrease the vessel numbers, thereby contributing to the white earlobe formation via paling the redness of skin in QYP, but its specific mechanism remains to be further clarified. In addition to TP63, CDH1 as an anti-angiogenesis gene, its was identified as a hub node gene by PPI, suggesting its potential pivotal roles in regulating angiogenesis thereby affecting the white/red earlobe formation. Collectively, our findings advance the existing understanding of the white earlobe formation, as well as provide new clues to understand the molecular mechanism of chicken white/red earlobe color formation.

This work was supported by grants from Science and Technology Planning Project of Guangzhou City (201604020007; 201504010017) and Guangdong Province (2014B020202012), and National Natural Scinence Foundation of China (31761143014). We also thank Endashaw Jebessa of Ethiopia for his kind edition on this manuscript.

Wei Luo performed the experiments, analyzed the data, wrote the paper, prepare tables and Fig. Jiguo Xu prepared the experiment samples, analyzed the data, and participated in the design of the study. Zhenhui Li, Haiping Xu and Jiaying Wang helped analyze the data. Shudai Lin and Hongjia Ouyang helped perform experiments. Qinghua Nie conceived of the study and participated in its design and coordination. Xiquan Zhang participated in the design of the study.

The authors declare to have no conflicts of interest.

W. Luo and J. Xu contributed equally to this work.