We used a local Chinese Baicheng You chicken population to screen candidate genes for feather colors. First, 30 adult chickens were selected from differentiated sub-populations selected for the three target feather colors (Figure 1), lavender, black, and yellow, with 10 samples for each sub-population (ADMIXTURE and PCA; Supplementary Figure S1). Blood samples were obtained from the brachial veins by standard venipuncture and were placed into centrifuge tubes containing an anticoagulating agent.

Genomic DNA was extracted from the 30 chickens using the standard phenol/chloroform method, and random interruptions were made using a Covaris ultrasonic crusher. We performed 150-bp paired-end re-sequencing using the Illumina HiSeq 2500/NovaSeq 6000 system according to the manufacturer’s protocols. We used fastp (v0.20.0) software to filter residual primers, adapters trimming, and lower-quality reads, and used the -q 20 -u 30 parameters (from the Fastp manual, q, --qualified_quality_phred indicate the quality value that a base is qualified, thus default 15 means phred quality ≥ Q15 is qualified) for the quality control of the 150-bp paired-end raw reads (Chen et al., 2018). If the quality value of the per-base or per-read is less than 20, this is considered to indicate low quality. If the percentage of low-quality bases is 30%, the reads with too many N bases are excluded. If the content of N bases is greater than n, the read/pair will be discarded, and the default value is 5, such as when one has a high quality score (>Q30) and the other has a low quality score (<Q15). To reduce false corrections, fastp only performs a correction if the total mismatch is below a given threshold T (T ¼ 5) (Chen et al., 2018). The clean reads were mapped to the chicken genome (https://www.ncbi.nlm.nih.gov/, assembly: GCA_016700215.2) using BWA-MEM (v0.7.17) with default parameters, except for -t 4 -R (-t is thread, -R is set the reads header, and it is split with \t) option (Li et al., 2008). We sorted the alignment bam files using the Samtools (v-1.11) software package, and duplicate reads were removed using the picard (http://broadinstitute.github.io/picard/) tools MarkDuplicates. Genome Analysis Toolkit (GATK,v-4.2) (https://gatk.broad institute.org/) was then used for downstream processing and variant calling. Using the HaplotypeCaller algorithm of the GATK pipeline in the genomic Variant Call Format (gVCF), we obtained the genotype likelihoods at each site in the reference genome for each individual. We called the SNPs and indels using the HaplotypeCaller algorithm. Individual gVCF files were then merged into a multi-individual VCF using the GenotypeGVCFs tool of the GATK pipeline. It is imperative to filter raw SNP candidates in the genotyping workflow because it allows the shrinking of false-positive calls due to biases in the sequencing data (Jiang et al., 2021; Xiangtao Liu et al., 2013). We called SNPs and indels using the SelectVariants algorithm. Next, these variants were used as input for hard filtering in the GATK pipeline based on six statistics to identify SNPs or indels: QUAL < 30.0; QualByDepth (QD) < 2.0; FisherStrand (FS) > 60.0; RMS MappingQuality (MQ) < 40.0; MappingQualityRankSunTest (MQRankSum) < -12.5; and ReadPosRankSumTest (ReadPosRankSum) < -8.0 (Li et al., 2008). ClusterWindowSize and clusterSize were set to 10 and 3, respectively. Generally, consecutive SNPs tend to provide a high proportion of false-positive results in terms of selection signatures associated with plumage color (Margot E. Bowen et al., 2011).

The cross-population composite likelihood ratio (XP-CLR) method is mainly based on the gene Site Frequency Spectrum (SFS) principle, which is based on the difference in multilocus allele frequency between two populations and is used in selection signal detection (Maria Ines Fariello 2013; Chen et al. 2010). These methods have been widely used to discover the loci and genes related to difference in traits, such as hair length in Tianzhu white yaks (Bao et al., 2022), body size in dogs (Lyu et al., 2021) and ponies (Asadollahpour et al., 2020), and meat quality in Angus cattle (Taye et al., 2018). We used XP-CLR (v1.1.2) (Chen et al., 2010) statistical methods to detect genetic differentiation in Baicheng You chicken plumage colors (Pavlidis et al., 2013). The values of XP-CLR were calculated using the python script XPCLR, downloaded from Github (https://github.com/hardingnj/xpclr). The corresponding parameters were set as follows: maximum SNPs = 600, ld values = 0.95, window size = 40,000 (Lee et al., 2018). The regions with XP-CLR values in the top 1% were considered candidate regions, and genes with overlapped candidate sweeps were considered candidate genes. Six feather color models were used to identify candidate direct selection regions: 1) black vs. lavender; 2) black vs. yellow; 3) yellow vs. lavender; 4) black vs. yellow + lavender; 5) lavender vs. black + yellow; and 6) yellow vs. black + lavender. In addition, the genes were annotated using the chicken genome (https://www.ncbi.nlm.nih.gov/, assembly: GCA_016700215.2) and NCBI databases. Variant effect predictor (VEP) was used to annotate gene variants according to their functional categorization, including the following categories: up- and downstream gene, frameshift, 3- and 5-prime UTR, synonymous, missense, start lost, intronic, and splice region. We used the R qqman package to create a Manhattan plot of XP-CLR for the association analysis.

After quality control, a total of ∼8.9 million high-quality paired-end reads were obtained from the 30 chickens (267.73 Gb clean base). The average sequencing depth per individual was 8.11×, and the average genome coverage was 98.62%.

A total of 18.37 million SNPs and 2.61 million indels were retained for the association analysis. In addition, where more than 3 SNPs were clustered within a 10-bp window, they were all considered false positives and deleted. In the subsequent analyses, we used only biallelic SNPs on the chromosomes.

The candidate genes for the three kinds of feather colors were identified using XP-CLR statistical methods in Baicheng You chicken. The total numbers of genes with positive selection signatures detected in the six groups were as follows: 934 genes (black vs. lavender), 3620 genes (black vs. yellow), 989 genes (yellow vs. lavender), 685 genes (black vs. yellow + lavender), 836 genes (lavender vs. black + yellow), and 820 genes (yellow vs. black + lavender), including cell growth, differentiation, migration, and immune regulation (Supplementary Table S2–S7). XP-CLR scores are presented in Figure 2.

First, to obtain the genes associated with black plumage, the scanning regions of top 1% of black vs. lavender, black vs. yellow, and black vs. yellow + lavender models were extracted, respectively (Figures 2A,B,D), and the intersection of the three sets of data was taken. We obtained the overlap windows at chromosome 13: 18940001–18980000 bp, which overlapped with EGR1 gene.

Second, we extracted top 1% of sweep regions of black vs. lavender, lavender vs. yellow, and lavender vs. black + yellow models, respectively (Figures 2A,C,E), and the same method was used to obtain the intersection of the three groups of data. One window was found at chromosome 7: 4800001–4840000 bp, harboring RAB17, and we also identified a window at chromosome 7: 4820001–4860000 bp, containing MLPH. Another window was located at chromosome 1: 65940001–65980000 bp, which overlapped with the SOX5 gene. We thus obtained MLPH, RAB17, and SOX5 as genes that could be associated with lavender plumage color.

We detected the top 1% of sweep regions of black vs. yellow, lavender vs. yellow, and yellow vs. black + lavender models, respectively (Figures 2B,C,F), and again took the intersection of three sets of data. The overlap window at chromosome 1: 18240001–18280000 bp, containing the GRM5 gene.

The lavender phenotype has been found to be related to MLPH and Rab17 (Roulin and Ducres, 2013). SOX5, as a transcription factor, indirectly regulates the synthesis and transport of melanocytes (Min et al., 2022; Stolt et al., 2008; Wang et al., 2016). EGR1 is expressed in mouse hair follicles (Mcmahon et al., 1990). GRM5 is a candidate gene in the Polverara chicken (black vs. white) (Mastrangelo et al., 2020). Thus, the results showed that MLPH, RAB17, SOX5, EGR1, and GRM5 genes were associated with lavender, black, and yellow plumage colors in Baicheng You chicken.

A total of 258 SNPs and 28 indels were identified in MLPH (Supplementary Table S8). In the coding regions of the MLPH gene, we detected eight synonymous, two non-synonymous nucleotide substitutions, and one deletion. Interestingly, we found an A99G substitution in exon 1 and one deletion (AG861A) in exon 6 that were present in lavender feather chickens.

RAB17 possessed a total 179 SNPs and 6 indels (Supplementary Table S9). Five synonymous and four non-synonymous nucleotide substitutions were detected in the coding region of RAB17 gene, and G514A (Val172Ile) substitution in exon 5 appeared in lavender plumage color (Table 1).Three synonymous and two non-synonymous nucleotide substitutions were detected in exon 2 of the EGR1 gene (Supplementary Table S10).

The remaining two plumage color loci at chromosome 1 exhibited the strongest association signals within the intronic regions of GRM5 and SOX5, and they are obvious candidates for controlling normal variation in chicken plumage colors (Supplementary Tables S11, S12). However, GRM5 lies 33.82 kb upstream of TYR, which is also involved in pigmentation, which could affect pigmentation by regulating the expression of TYR (Sandra Beleza et al., 2013; Smyth et al., 2006). Therefore, we hypothesized that the non-coding DNA variation of GRM5 and SOX5 might affect pigmentation by regulating the expression of neighboring genes or influence itself expression by regulating the enhancer activity of the genes themselves. Candidate genes based on the SNPs and amino acid substitutions in black, lavender, and yellow feathers in Baicheng You chicken are summarized in Table 1.

In our study of the three feather colors of Baicheng You chicken, we focused specifically on the loci that may play a more important role in the evolution of feather color through natural, sexual, and artificial selection. Here, we used comparative genome analysis between black, yellow, and lavender sub-populations to identify the genetic basis underlying the variation in feather color among Baicheng You chicken. In this study, several previously reported genes were found to be involved in lavender feather-related traits, such as MLPH, RAB17, and SOX5. In previously published articles, non-synonymous mutation (c. 103 C>T) (Vaez et al., 2008) and (c. 1909 A>G) (Xu et al., 2016) of MLPH resulted in the development of the lavender or light-gray plumage color phenotype in chicken. In this study, we did not find SNP (c.1909 A>G), but a frameshift deletion (c. 861 indel G) was located at exon 6 of MLPH, which leads to a completely different protein. Vaez et al. found that the absence of exon 6 in MLPH gene transcript in lavender phenotype chickens (Vaez et al., 2008), and frameshift deletion in exon 5 of MLPH gene has similarly been reported in rabbit (Oryctolagus cuniculus) breeds (Fontanesi et al., 2014). Deletions have also been found in leaden mice and cat (Ishida et al., 2006). The variations in MLPH gene could cause transport defect melanosomes in melanocytes, resulting in the dilution of coat or plumage colors in cats (Ishida et al., 2006), dogs (Bauer et al., 2018; Drögemüller et al., 2007; Van Buren et al., 2020), cattle (Li et al., 2016), sheep (Posbergh et al., 2020), rabbits (Demars et al., 2018; Jia et al., 2021), and minks (Cirera et al., 2013), as well as Griscelli syndrome type-3 in humans (Al-Mousa et al., 2016; Çağdaş et al., 2012; Gironi et al., 2019; Jo et al., 2020; Ménasché et al., 2005). McMurtrie et al. reported that RAB17 is an attractive candidate for leaden mice, though failed to identify any RAB17-coding region mutations, because it was expressed in polarized epithelial cells (E. B. McMurtrie et al., 1997). A non-synonymous substitution (G514A) in exon 5 of RAB17 gene by whole genome re-sequencing was detected, resulting in change in the amino acid Val172Ile in this study.

The transcription factor SOX5 has been detected in the early migrating neural crest of chicken (Perez-Alcala et al., 2004). There are two explanations for the regulation of transcription factor SOX5 on melanocytes. A 7.6-kb non-coding deletion near SOX10 was found to be associated with pale yellow feathers (Zhu et al., 2022), and SOX5 played a role in melanocyte lineage by regulating the activity of transcription factor SOX10 (Min et al., 2022). The other was up-regulation of the expression of TYR through MITF, thereby reducing the synthesis of melanin (Wang et al., 2016).

In chicken, MC1R is a single exon gene without intron on chromosome 11, and it is associated with the traditional feather color Extension locus (E). Although MC1R is among the genes that have been extensively studied for their association with feather color, we did not identify it except for in the black vs. yellow + lavender model. We had the major gene EGR1 of black plumage by XP-CLR algorithm, and we hypothesized that the variation loci T1297C (Ser433Phr) in EGR1 could be associated with the synthesis of eumelanin. EGR1, a member of the immediate early family, is an important nuclear transcription factor (Gómez-Martín et al., 2010). EGR1 mRNA was highly expressed in hair follicles in the process of mouse embryogenesis (Mcmahon et al., 1990). EGR1 could also up-regulate the expression of TYR by binding MC1R, leading to the synthesis of large amounts of eumelanin.

In this study, GRM5 was associated with yellow plumage color. Our results coincided with those of Mastrangelo et al. (2020). In recent years, studies of human skin and eye color have found that the skin color loci at 11q14.3 is strongly correlated with the intronic regions of GRM5, and non-synonymous mutation in TYR (S192Y) is in linkage disequilibrium with the most significant SNPs in GRM5 (Sandra Beleza et al., 2013). Kaustubh Adhikari et al. have also demonstrated the presence of multiple independent signals of association in GRM5/TYR (Adhikari et al., 2019; Miao et al., 2017). GRM5 may affect pigmentation by regulating the expression of TYR.