A quantitative trait locus on chromosome 2 was identified that accounts for a substantial proportion of phenotypic variance of the yellow plumage color in chicken.

Chicken plumage color is an important economical trait in poultry breeding, as triple-yellow indigenous broilers are preferred over western commercial broilers in the Chinese market. However, the studies on the pigmentation of plumage coloration are relatively rare at present. Here, we performed a genome-wide mapping study on an F2 intercross, whose 2 founders were one hybrid commercial line "High Quality chicken Line A" that originated from the Anak red chicken and one indigenous line "Huiyang Beard" chicken that is a classical "triple-yellow" Chinese indigenous breed. Moreover, we used an automatic colorimeter that can quantitatively assess the colorations in L∗, a∗, and b∗ values. One major quantitative trait locus (QTL) on chromosome 2 was thus identified by both genome-wide association and linkage analyses, which could explain 10 to 20% of the total phenotypic variance of the b∗ measurements of the back plumage color. Using linkage analysis, 2 additional QTL on chromosome 1 and 20 were also found to be significantly associated with the plumage coloration in this cross. With additional samples from Anak red and Huiyang Beard chickens as well as pooled resequencing data from the 2 founders of this cross, we then further narrowed down the QTL regions and identified several candidate genes, such as CABLES1, CHST11, BCL2L1, and CHD22. As the effects of QTL found in this study were substantial, quantitatively measuring the coloration rather than the descriptive measurements provides stronger statistical power for the analyses. In addition, this major QTL on chromosome 2 that was associated with feather pigmentation at the genome-wide level will facilitate the future chicken breeding for yellow plumage color. In conclusions, we mapped 3 associated QTL on chromosome 1, 2, and 20. The candidate genes identified in this study shed light in the genetic basis of yellow plumage color in chicken.

Introduction

In China, chicken plumage color is an important economic trait, as most native breeds have characteristic “triple-yellow” appearance (yellow plumage, beak, and claw) and they are favored by Chinese consumers. Although colorations of plumage were used to be categorized descriptively, its underlying genetic basis is complex, which included the biosynthesis of eumelanin and pheomelanin in melanocytes (Emaresi et al., 2013).

In recent years, genome-wide association studies (GWAS) were widely applied in coat coloration of animals, including Iranian Markhoz goats (Nazari-Ghadikolaei et al., 2018), mice (Miao et al., 2017), horses (Kim et al., 2017), and dogs (Hedan et al., 2019), to screen novel candidate regions. For the black and brown colors of chicken plumage, previous studies adopted GWAS to discover significant SNP on chromosomes 1, 2, 3, and 5 (Yang et al., 2017) and on chromosomes 1, 2, 3, 4, 6, 8, 19, 24, and 33 (Park et al., 2013). For the chocolate color in chicken plumage, recent study using linage analysis (LA) reported that a missense mutation in TYRP1 alters the melanosome structure, which caused the pigment of chocolate color (Li et al., 2019).

So far, many well-known genes including ASIP, TYR, MC1R, TYRP1, SLC24A4, MITF, and KIT have been identified to play vital roles in melanocyte differentiation, proliferation, migration, survival, and morphology (Sulem et al., 2007, Duffy et al., 2010, Minvielle et al., 2010, Roulin and Ducrest, 2013). In feather follicles, the binding of melanin-stimulating hormone α-MSH to Mc1R could increase the intracellular cAMP level, and this promotes MITF expression and further increases the activity of eumelanin-related enzymes such as Tyr and Tyrp1 (Garcia-Borron et al., 2005). Besides, the expression of ASIP in the Japanese quails was higher in the yellow stripes of their back skin than that in the black stripes (Inaba et al., 2019). The binding of ASIP to Mc1R could induce the production of pheomelanin (Sakai et al., 1997, Walker and Gunn, 2010).

Owing to the studies on the pigmentation of pheomelanin (yellow, red) being relatively rare at present, we attempted to quantify the coloration of chicken back plumage. Based on an F2 intercross population, we mapped 3 associated quantitative trait loci (QTL) on chromosome 1, 2, and 20. Then, using genotypes from additional birds and whole genome pooled resequencing data of the founder lines, QTLs were further narrowed down and candidate genes were identified. Our results provide new insights into the genetic basis of yellow plumage color in the chicken.

Materials and methods

Experimental Population

The F2 intercross population and the corresponding genotyping data have been analyzed for its growth traits in the study by Sheng et al. (2013). Detailed descriptions can be found in that study. In brief, the founders were one hybrid commercial line “High Quality chicken Line A” (HQLA) that originated from the Anak red chicken and one indigenous line “Huiyang Beard” (HB) chicken that is a classical “triple-yellow” chicken breed. By reciprocal crossing of the founder lines for 2 generations, a total of 800 F2 birds were produced. Additional samples from Anak red and HB breeds (100 each) were collected from the National Chickens Genetic Resources (Jiangsu Province) (http://www.genebank.org.cn/), which is an off-site conservation base for 29 selected breeds.

Phenotyping and Genotyping

In this study, an automatic colorimeter TCP2 (Beijing Oyke Opto-Electronic Instrument Co., Ltd.) was used to quantitatively assess the back plumage color. The color is quantified using values a∗ (red to green axis; ABPC) and b∗ (yellow to blue axis; BBPC). Light to dark (LBPC) is determined as L∗ on a separate numerical scale (from black = 0 to white = 100). The measurements were performed according to Watson et al. (2018). Birds of the F2 generation were c`onsecutively tested every other week from 8 to 12 wk of age. Four different points on the back plumage of each chicken were measured to obtain an average value.

Blood samples were collected from a total of 585 birds from the F2 cross and 200 birds from the Anak red and HB chickens. Genotyping for all samples were performed using Illumina Chicken 60K SNP Beadchip (Groenen et al., 2011). Quality controls were performed with the following standards: individual data were excluded when their call rates were below 0.9; single-nucleotide polymorphisms (SNP) were removed as a result of call rates < 0.9, minor allele frequency < 0.05, or undetermined positions on the chromosome. SNP on Z chromosome were kept only when their genotypes were homozygous in hens. A final of 44,177 SNPs from the 494 birds of F2 cross were used for the next-step analyses.

Whole-Genome Linkage Analysis

Genetic map was constructed by using CRI-MAP, which was about 3,068 cM in total. Linkage mapping was carried out with least squares regression as described by Haley et al., 1994. With customized R scripts, a multiple linear regression model was used with “sex” and “batch” as the fixed effects. Effects of putative QTL were estimated at each 1-cM interval across the genome, and F statistics were calculated accordingly. One percent and 5% genome-wide thresholds were determined by 1,000 times permutation tests as described by Churchill and Doerge, (1994). To further explore the independent QTL, the most significant QTL would be added to the model as covariate(s) to perform another round of analysis, until no significant result was obtained.

The proportion of the residual phenotypic variance that was explained by the detected QTL was calculated by the following equation:

where, is the residual mean square of the reduced model (i.e., model that included “sex” and “batch” effects, also other QTL as cofactors), MSF is the residual mean square of the full model (i.e., model that included all covariates, other QTL as well as the additive and dominance effects of the tested position), and MSR is the residual mean square of the reduced model (model that only fit all fixed effects but not any QTL for single-QTL analysis). The methods used for linkage analysis in this study were the same as the study on growth traits of the same population; details can be found in the study by Sheng et al. (2013).

Genome-Wide Association Analysis

A mixed linear model was used to perform association analysis in this study. While “sex” and “batch” were included as the fixed effects, the polygenic effect was also contained in the model as a random effect to correct for population stratifications. The association analyses in this study were all based on this model and were implemented by the R package – GenABEL (Aulchenko et al., 2007). A two-step score test was performed. By using the command “ibs,” the kinship matrix of the experimental population was constructed. Then, the matrix was implemented into the function “polygenic” to obtain the phenotypic residuals that had been accounted for 2 fixed effects and one polygenic (random) effect. The narrow sense heritability was also a direct output of the function “polygenic.” Finally, the obtained residuals were used as dependent traits in a simple regression for single-marker or haplotype-based association tests.

The proportion of the phenotypic variance explained by each of the significant SNP was calculated as follows:where, p and q are the frequencies of the 2 alleles at the tested locus and Vary is the variance of the tested trait. As a set of significantly associated SNP will be identified for the feather pigmentation at different weeks of age, we only estimated the phenotypic variances explained by the most significant SNP. Thus, is the marker genotype effect of the most significant SNP to the trait.

Phasing and Haplotyping

The selected 3 peak SNP from 3 populations were first phased using fastPHASE 1.2, separately (Scheet and Stephens, 2006), and then, the frequencies were calculated. Haplotypes that had frequencies more than 5% were retained for further analysis. Using a customized R script, the 2 core haplotypes identified were plotted against their corresponding residual phenotypes, which were obtained in the reduced model that only accounted for 2 fixed effects.

FST Calculation Using Whole-Genome Resequencing Data

Whole-genome resequencing on pooled samples from the 2 founder lines (15 HB and 16 HQLA birds) were performed. Four hundred base pair paired-end sequencing libraries were prepared for the next whole-genome resequencing using the Illumina HiSeq 2,000 platform (BGI, Shenzhen). Clean raw reads provided by BGI were mapped onto the chicken reference genome (Gallus_gallus-5.0) using the command “mem” in BWA (Houtgast et al., 2018) and were then sorted by SAMtools (Li et al., 2009); the resulting Binary Alignment/Map files were next processed by Picard 2.8.0 to remove the PCR duplicates; a further recalibration on the base quality score of the aligned data before SNP calling was performed using “BaseRecalibrator” in GATK (Ni and Stoneking, 2016); SAMtools again was used to generate mpileup files for PoPoolation2 to calculate FST over 1 kb sliding window with 50% overlaps. The settings used for PoPoolation2 (Kofler et al., 2011) were minimum count 3, minimum coverage 10, maximum coverage 100, and minimum coverage fraction 1 (Lillie et al., 2017). Threshold was then set to 0.43 as this was the 1% top FST value across the genome in this study. The r2 values of linkage disequilibrium decreased over shorter distances (an average of 3.5 kb) in HB chicken than in Anak red (Peng, 2015). Therefore, we clustered significant regions (FST > 0.43) that were no more than 3.5 kb apart.

Selection of Candidate Genes in the QTL

For GWAS results, the SNP location was referred to the Gallus_gallus-5.0 reference genome, and the genes within the region that 0.35 Mb upstream and downstream of the associated SNP were considered as the candidates. For LA results, genes in the peak regions were considered as candidates.

Results

Phenotypic Summary

The statistical summaries of the L∗, a∗, and b∗ values of back plumage color through 8 to 12 wk of age (ABPC8-12, BBPC8-12, and LBPC8-12) were shown in Supplementary Table 1. The coefficient of variation (CV) values were all at relatively low levels (0.13–0.24), which suggested our measurements were rather stable within this population. Based on the measurements of ABPCs and BBPCs, the back plumage colors were of different degrees of red and yellow, which were in consistent with the fact that they were hybrids of Anak (red) and Huiyang Beard (yellow) chickens.

Pearson's correlation coefficients were calculated between each pair of the recorded phenotypes by “cor.test” in R. ABPC were negatively correlated with BBPC and LBPC (Table 1). The estimated r values between LBPC and ABPC or BBPC were at intermediate to high levels (0.4–0.8), whereas those between ABPC and BBPC were rather low. These are in accordance with the P values; the lower the r values, the higher the P values. In summary, both ABPC and BBPC were highly significantly and closely correlated with LBPC. However, little or low correlations were observed between ABPC and BBPC.

Table 1

Correlations between pairs of traits recorded in this study.

	ABPC8	ABPC10	ABPC12	BBPC8	BBPC10	BBPC12	LBPC8	LBPC10	LBPC12
ABPC8	1	0.50	0.36	−0.03	−0.04	−0.16	−0.52	−0.44	−0.35
ABPC10	9.69E-32	1	0.71	−0.29	−0.13	−0.22	−0.57	−0.75	−0.63
ABPC12	9.65E-17	2.02E-76	1	−0.27	−0.28	−0.28	−0.48	−0.66	−0.81
BBPC8	5.13E-01	3.21E-11	7.99E-10	1	0.64	0.59	0.66	0.58	0.52
BBPC10	3.72E-01	3.92E-03	1.73E-10	5.26E-58	1	0.76	0.53	0.61	0.58
BBPC12	3.19E-04	1.22E-06	3.93E-10	1.50E-47	1.33E-91	1	0.54	0.63	0.70
LBPC8	4.12E-36	1.57E-42	1.04E-29	5.22E-63	1.52E-36	1.53E-38	1	0.76	0.67
LBPC10	2.87E-24	1.03E-90	5.57E-62	7.69E-46	1.19E-50	2.36E-55	9.40E-92	1	0.85
LBPC12	2.78E-15	3.49E-55	9.08E-114	1.82E-35	1.48E-45	1.87E-72	3.79E-65	5.34E-139	1

Diagonal and upper triangle: Phenotypic correlations. Lower triangle: P values of the phenotypic correlations.

Three Independent QTL Were Found to be Significantly Associated With the L∗ and b∗ Values of Back Plumage Color in This Population

Given the recorded pedigree of this population, genetic map was constructed and both GWAS and LA were performed (all results were listed in Supplementary Table 1 significant results were summarized in Table 2). For the ABPC measurements, only a single significant SNP was found that was associated with ABPC10 in GWAS. However, based on the increasing order of the corresponding P values, the following SNP were identified on chromosome 20. In addition, this locus was not significantly associated with ABPC10 using LA, and the locus with the highest F value in LA was on the other chromosome – GGA10 (Supplementary Table 1). Therefore, no significant results were found for ABPC measurements in this study.

Table 2

Loci that were significantly associated with the plumage color on birds' back in this population.

Traits (h2)	QTL	Chromosome	Peak positions	Test statistics (F/P values)	Peak regions	Variance explained	Covariates	Methods
LBPC8 (43.7%)	QTL2	2	103.2 Mb (rs14228181)	1.25E-09∗∗	103.2-106.7 Mb	7.85%	/	GWAS
	QTL2	2	102.4 Mb (222 cM)	18.56∗∗	88.9-108.7 Mb (201-235 cM)	6.73%	/	LA
LBPC10 (52.9%)	QTL2	2	103.2 Mb (rs14228181)	2.91E-08∗∗	103.2-106.7 Mb	6.65%	/	GWAS
	QTL2	2	101.7 Mb (221 cM)	24.63∗∗	85.5-111.9 Mb (199-240 cM)	8.83%	/	LA
	QTL1	1	54.1 Mb (139 cM)	13.64∗∗	53.6-56.1 Mb	4.51%	chr2:221 cM	LA
LBPC12 (46.2%)	QTL2	2	106.7 Mb (rs317554548)	3.19E-07∗	/	5.59%	/	GWAS
	QTL2	2	101.1 Mb (220 cM)	17.95∗∗	91.3-108.1 Mb (207-233 cM)	6.48%	/	LA
	QTL5	5	53.9 Mb (149 cM)	11.04∗∗	/	3.70%	chr2:220 cM	LA
BBPC8 (44.0%)	QTL2	2	103.2 Mb (rs14228181)	8.06E-13∗∗	93.8-106.7 Mb	12.08%	/	GWAS
	QTL2	2	101.7 Mb (221 cM)	36.91∗∗	83.1-113.2 Mb (195-242 cM)	12.85%	/	LA
	QTL20	20	10.2 (63 cM)	11.96∗∗	9.9-10.4  Mb (61-65 cM)	3.77%	chr2:221 cM	LA
BBPC10 (56.7%)	QTL2	2	103.2 Mb (rs14228181)	1.27E-11∗∗	93.8-107.6 Mb	11.50%	/	GWAS
	QTL2	2	102.9 Mb (223 cM)	57.86∗∗	76.1-116.4 Mb (189-249 cM)	19.02%	/	LA
	QTL20	20	10.3 Mb (64 cM)	20.01∗∗	9.3-13.3 Mb (56-75 cM)	5.92%	chr2:223 cM	LA
	QTL1	1	50.5 Mb (126 cM)	13/31∗∗	49.5-53  Mb (123-135 cM)	3.66%	chr2:223 cM; chr20:64 cM	LA
BBPC12 (55.7%)	QTL2	2	103.2 Mb (rs14228181)	4.61E-11∗∗	100.1-108.3	10.28%	/	GWAS
	QTL2	2	102.4 Mb (222 cM)	45.62∗∗	81.3-112.3 Mb (194-241 cM)	15.54%	/	LA
	QTL20	20	10.8 Mb (67 cM)	10.98∗	10.6-10.9 Mb (66-67 cM)	3.35%	chr2:222 cM	LA

∗Denotes 5% genome-wide significance, ∗∗denotes 1% genome-wide significance.

Abbreviations: GWAS, genome-wide association studies; LA, linage analysis; QTL, quantitative trait locus.

Using GWAS, 6 QTLs were identified, and each was significantly associated with one measurement of BBPC and LBPC. Twelve QTLs were found to be associated with the 6 traits of BBPC and LBPC by LA. As shown in Table 2, these detected QTLs clustered to 4 locations on 4 different chromosomes in the genome. We then named each QTL by combining “QTL” and the number of chromosome where it was located.

One major QTL on chromosome 2 (QTL2) were not only identified by both methods but were also associated with LBPC and BBPC across 8 to 12 wk of age. For LBPC, QTL2 explained a large proportion of phenotypic variances (5.5–8.8%). An even larger proportion of phenotypic variances of BBPC were explained by QTL2 (10–19%; Table 2, Figures 1A and 1C). After correcting for the effect(s) of the most significant SNP as covariate(s) in the model, additional 3 QTLs on chromosomes 1, 5, and 20 were found by LA, whereas no significant results were detected by GWAS. As QTL5 was identified to be associated with LBPC12 alone and it spanned only one centimorgan, we excluded it as an unreliable result from this study. Therefore, another 2 independent QTL20 and QTL1 were found to be associated with BBPC by LA (Figures 1B and 1D).

Figure 1

Genome-wide association and linkage analysis for LBPC10 and BBPC10 traits. (A) Manhattan plot for the genome-wide association study for LBPC10 traits. (B) The identification of QTL related with LBPC10 trait by linkage analysis. (C) Manhattan plot for the genome-wide association study for BBPC10 traits. (D) The identification of QTLs related with BBPC10 trait by linkage analysis. The Manhattan plot indicates the −log10 (observed P values) for genome-wide SNP (y-axis) plotted against the irrespective positions on each chromosome (x-axis).

Identifying the Core Haplotypes Underlying the Major QTL on Chromosome 2

Compared with the estimated heritability of each BBPC traits, the peak SNP of QTL2 accounted for about 20% of the additive genetic variances, which suggested QTL2 was a major QTL that was associated with measurements of BBPC in this population. We then further define the haplotypes using the top 3 peak SNPs found by GWAS for trait LBPC10. Two main haplotypes “AGG” and “GAA” were identified, whose frequencies took up nearly 80% of all the haplotypes found. By using residual phenotypes that were obtained by corrected for sex and batch effects, “AGG” homozygous individuals had the highest b values, whereas “GAA” homozygous individuals had the lowest b values (Figure 2). We further examined the frequency distributions of “AGG” and “GAA” in additional 200 Anak red and HB birds (100 birds each; Table 3). “AGG” that increased the extent of yellowness was not found in the Anak chicken, whereas the frequency of “GAA” was substantial (0.78). Moreover, “GAA” took up significantly smaller proportion in the HB chicken.

Figure 2

The distribution of LBPC10 measurements on the different core haplotypes. “AGG” and “GAA” were the 2 main haplotypes. Letters in yellow mean the corresponding haplotypes have an effect to increase L∗ values, while letters in blue denote their effects in decreasing the L∗ values.

Table 3

Frequency distributions of core haplotypes from QTL2 in 3 populations.

Populations	AGG	GAA	others1
F2 population	0.41	0.36	0.23 (3)
Anak red	0	0.78	0.22 (2)
HB	0.07	0.24	0.69 (4)

Numbers of haplotypes that had frequencies more than 5% were shown in brackets.

QTL Were Further Narrowed Down Using Whole-Genome Pooled Resequencing Data of the Founder Lines

The candidate genes listed in Supplementary Table 1 were determined based on the position of the associated SNP identified by GWAS and the promising regions identified by LA. Whole-genome FST values were obtained by comparing the HB and HQLA pooled sequencing data, and the FST profiles of the QTL regions were shown in Figure 3. Almost 200 significant FST clusters were defined within the 3 identified QTLs (Supplementary Table 1). On QTL1, the candidate genes including MIEF1, MGAT3, TNRC6B, and CHST11 were overlapped with FST clusters. The CHST11 is a key enzyme in the biosynthesis of chondroitin sulfates. On the core haplotypes of the QTL2, the genes that contain the associated SNP were CABLES1, TAF4B, GREB1L, ROCK1, THOC1, COLEC12, and ANKRD29, besides, 4 candidate genes including ZNF521, LAMA3, TMEM241, and CABLES1 were overlapped with the significant FST clusters. Furthermore, there were 104 candidate genes within the region of the QTL20. Of these candidates, BCL2L1 and CDH22 were overlapped with the significant FST clusters. Combined with the implication of the relevant published references, the promising genes were emphasized and summarized in Table 4.

Figure 3

The distribution of the whole-genome FST clusters within the QTL. The arrows represent FST clusters in the core haplotypes. 1% top FST clusters across the genome was defined as the threshold set to 0.43.

Table 4

The functions of candidate genes identified in this study.

QTLs	Gene	Position	Function	Reference
1	CHST111	54433922-54593321	Catalyze the transfer of sulfate to position 4 of the GalNAc residue of chondroitin.	Zhou et al., 2016
2	DSG4	106052985-106078954	Haplotypes related with the color of the goat coat.	E et al., 2016
2	CABLES11	102881201-102953764	An adapter protein that links cyclin-dependent kinase with nonreceptor tyrosine kinases and regulates the activity of Cdks by enhancing their Y15 phosphorylation.	Huang et al., 2017
2	GATA6	102636810-102654409	Promotes hair follicle progenitor cell renewal.	Wang et al., 2017
2	GREB1L	102280644-102409696	A coactivator for retinoic acid receptors.	Brophy et al., 2017
2	YES1	101881440-101930073	Mediates the biological process of tyrosine.	Patel et al., 2013, Ta et al., 2018
20	EDN3	10799441-10816861	Variation in EDN3 resulted in dermal hyperpigmentation phenotype or fibromelanosis	Dorshorst et al., 2011, Lillie et al., 2017, Shinomiya et al., 2012
20	BCL2L11	10015061-10033228	Homologous to BCL2 gene related with gray hair in mouse	Seiberg et al., 2013
20	CDH221	10647698-10705460	The expression was associated with melanoma metastasis	Piche et al., 2011

Genes overlapped with the significant FST clusters.

Discussion

Over the last few decades, researchers attempted to develop different methodologies and approaches to measure the coloration of animal appearance, such as Munsell tables (Peiponen, 1992), photographic analysis (Fitze and Richner, 2002, Figuerola and Senar, 2010), digital colorimeters (Figuerola et al., 1999, Juan Carlos et al., 2003), or reflection spectrometers (Bennett et al., 1996, Hunt et al., 1997, Alonsoalvarez et al., 2004). In our study, we quantified the plumage colorations using an automatic colorimeter. The measurements of plumage coloration can be conducted in 2 ways: one way was measuring the coloration directly on the bird (Siitari and Huhta, 2002, Eaton and Lanyon, 2003, Alonsoalvarez et al., 2004) and the other way was collecting feathers for measurement in the laboratory (Zuk and Decruyenaere, 1994, Bennett et al., 1997, Keyser and Hill, 1999, Siefferman and Hill, 2003). Earlier studies have demonstrated that the measured values obtained directly on the animal or from collected feathers, were nicely correlated (Quesada and Senar, 2006).

Traditionally, the plumage color is regarded as a qualitative trait; in this study, we developed the measuring method to obtain quantified values of the plumage color for GWAS and LA, and the results showed the QTL distributing on chromosome 1 (49.5–56.1 Mb), 2 (76.1–116.4 Mb), and 20 (9.3–13.3 Mb). Besides, it is noticeable that the proportion of phenotypic variance explained by QTL2 was substantial, which we speculated was attributed to the improvement of the measurement method.

Our results indicated that the major genes were likely located in QTL2 and very likely in the core haplotypes. The core haplotypes constructed by using the 3 most significant SNPs (rs14228181, rs14227093, and rs14230934) were ranging from 102.3 Mb to 106.4 Mb where many candidate genes were located. Previous studies demonstrated that DSG4 was related with the goat white and black color (E et al., 2016), and GATA6 promotes hair follicle progenitor cell renewal (Wang et al., 2017). Besides, the dietary intake of tyrosine could influence the color expression of the hair coat of adult dogs (Watson et al., 2018), and the YES1 as a candidate of interested is a member of the Src family tyrosine kinases. Network results from geneMANIA with human database showed YES1 have the physical interaction with KIT. The KIT gene is confirmed to be responsible for several common colors in pigs and horses, such as Dominant White and spotting (Marklund et al., 1998, Fontanesi et al., 2010, Durig et al., 2017, Hoban et al., 2018). In addition, we found 2 associated SNPs were located in CABLES1 that was also overlapped with the significant FST clusters. It was reported that CABLES1 was a critical regulator of corticotrope proliferation (Roussel-Gervais et al., 2016), linked cyclin-dependent kinase with nonreceptor tyrosine kinases (Huang et al., 2017). In emydid turtles or male painted turtles, carotenoids promoted lutein access to increased yellow and red chroma, as well as reduced ultraviolet chroma and brightness (Steffen et al., 2019). The carotenoids in vivo can be transformed to retinoic acid by the metabolic reaction (Weaver et al., 2018). In this study, we found one associated SNP (rs315171337) located in the intron of the GREB1L, and this gene implicated in activation of retinoic acid receptor targets (Brophy et al., 2017).

Considering that the feather pigmentation possibly effected by multiple genes, the candidates on QTL1 and QTL20 require attention as well. As the candidates on QTL1, CACNA1I belong to the family members of voltage-gated calcium channels and involved in the potential role of neuronal calcium signaling in schizophrenia (Xie et al., 2018). CHST11 can catalyze the transfer of sulfate to position 4 of the GalNAc residue of chondroitin (Zhou et al., 2016). As the candidates on QTL20, EDN3 affects melanocyte proliferation and differentiation (Baynash et al., 1994). Studies on chickens have illustrated that the structural variation of EDN3 resulted in dermal hyperpigmentation phenotype or fibromelanosis (Dorshorst et al., 2011, Shinomiya et al., 2012, Li et al., 2017). Besides, BCL2L1 and CDH22 were overlapped with the significant FST clusters. BCL2L1 is homologous to BCL2, which was incorporated in the Color Genes database (http://www.espcr.org/micemut/) for gray hair in mice (Seiberg, 2013), and CDH22 expression was associated with melanoma metastasis (Piche et al., 2011).

SLC24A5 is a putative cation transporter belonging to the potassium-dependent sodium–calcium ion exchanger family. Association of the SLC24A5 with pigmentation has been established in zebrafish, mice, humans, and chickens (Lamason et al., 2005, Vogel et al., 2008, Liu et al., 2011).This gene was downregulated in yellow quail and negatively regulated by ASIP (Nadeau et al., 2008). In this study, we found many of its homologous genes in the candidates, including SLC41A2 (chr1), SLC35C2 (chr20), SLC52A3 (chr20), SLC12A5 (chr20), and SLC38A6 (chr5). However, their functions on pigment synthesis were rarely reported.

In summary, no significant QTL was found for ABPC measurements. For the BBPC and LBPC measurements, QTL2 was associated with LBPC and BBPC across 8 to 12 wk of age; QTL20 and QTL1 were 2 additional independent QTLs that were found to be associated with BBPC. Moreover, we constructed the core haplotypes using the 3 most significant SNPs, which narrowed the region to 102.3–106.4 Mb. Combined with the whole-genome FST values, we further narrowed down the QTL regions and found several candidate genes. Further analyses and experiments are required to validate their functions in the feather pigmentation in chickens.