1 INTRODUCTION

Pigmentation variation among birds is astonishing and is a consequence of the ability to create complex patterns through varied pigment deposition across the body as well as on individual feathers. Due to the usual homogeneous plumage appearance within species, it is challenging to pinpoint the underlying genetic cause for the observed phenotype diversity among species. The domestic chicken is a widely used model for genetic studies of the bewildering pigment variation in birds and has been used successfully so in the past (Andersson et al., 2020). Chickens exhibit a variety of within-feather pigmentation patterns of which five have been described to be the consequence of the interaction between a proposed major locus called Patterning (Pg) and other loci (Figure 1; Figure S1) (Smyth Jr, 1990). Two of these interacting loci, Extension (E) and Dark brown (Db), have been identified at the molecular level. E corresponds to the melanocortin 1 receptor gene (MC1R), one of the major pigmentation genes in vertebrates (see below). Db is caused by a 8.3 kb deletion 14 kb upstream of exon 1 in SOX10 on chromosome 1 (Gunnarsson et al., 2011). Pg is assumed to affect the distribution of dark pigment on individual feathers and to underlie penciling, which is characterized by elliptic bands of eumelanin on the feather (Figure 1b). Pg has been described as an incompletely dominant, autosomal gene located on chromosome 1 in close linkage with Db and Melanotic (Ml) (Carefoot, 1999; Moore & Smyth, 1972).

Phenotypic characterization. (a) Pigmentation phenotypes of the breeds used to generate the intercross, Light Brown Leghorn (left) showing wild-type color and Fayoumi (right) showing autosomal barring. Photograph credit: Laurence Verrier for the LBL male and David Gourichon for the Fayoumi female. (b) Within-feather pigmentation pattern in chickens. The proposed major Patterning (Pg) locus interacts with several other loci causing complex variation in melanin distribution across individual feathers. Wild-type alleles are indicated by “+.” E/MC1R—Extension/Melanocortin-1 receptor, Ml—Melanotic, Co—Columbian, Db—Dark Brown/SOX10. Modified after Andersson et al. (2020)

There are two distinct types of barring patterns in chicken, autosomal barring, and sex-linked barring. While sex-linked barring is adding a white bar on a pigmented background, autosomal barring is adding a black bar on a gold or silver background on individual feathers (Smyth Jr, 1990). Sex-linked barring is caused by the combined effect of regulatory mutation(s) and missense mutations in the tumor suppressor gene CDKN2A (Hellström et al., 2010; Schwochow Thalmann et al., 2017). Autosomal barring on the other hand has been suggested to be caused by the combined effect of the three autosomal loci Db, E, and Pg (Andersson et al., 2020; Carefoot, 1984,1999; Moore & Smyth, 1972) and is the trademark of some chicken breeds such as Fayoumi, Campine, and Westfälische Totleger.

MC1R is a key locus in pigmentation biology, and MC1R mutations are underlying pigmentation polymorphisms in many vertebrates. It codes for a G-protein-coupled receptor located in the plasma membrane of melanocytes. Upon activation by its agonist α-melanocyte-stimulating hormone (αMSH), MC1R goes through a conformation change triggering a signaling cascade, which eventually leads to an increase of intracellular cAMP levels, which in turn activates transcription of genes involved in pigment production. MC1R signaling promotes production of black/brown pigment (eumelanin), whereas low or no activity leads to the default production of red pigment (pheomelanin) (Garcia-Borron et al., 2005). In birds, associations of MC1R variants with pigmentation traits have been found both in domesticated species (Davila et al., 2014; Kerje et al., 2003; Nadeau et al., 2006; Ran et al., 2016; Takeuchi et al., 1996) as well as in wild birds (Baiao & Parker, 2012; Johnson et al., 2012; Lamichhaney et al., 2016; Mundy et al., 2004; San-Jose et al., 2015; Theron et al., 2001). There are at least six phenotypically distinct Extension/MC1R alleles described in chicken (Andersson et al., 2020), which are here listed according to dominance and decreasing amount of eumelanin: Extended black (E*E), Birchen (E*R), wild-type (E*N), brown (E*B), buttercup (E*BC), and wheaten (E*WH or E*Y depending on dominance). These alleles determine the distribution of black eumelanin across the body and MC1R is therefore considered to be involved in primary pattern formation (Andersson et al., 2020; Smyth Jr, 1990). Birds carrying the dominant E and R alleles usually exhibit very dark to black adult plumage while the remaining alleles create patterns of brown, salmon and wheat color in different body regions. In adult males, only two alleles (E and R) can be distinguished from wild-type, while the other variant alleles only cause distinct phenotypic effects in female plumage (Andersson et al., 2020; Smyth Jr, 1990).

In the present study, we generated an informative backcross involving the Fayoumi breed showing autosomal barring (expected genotype Pg/Pg, Db/Db, and R/R at the E/MC1R locus) and an inbred line of Light Brown Leghorn (LBL) showing no feather pattern and assumed to be wild-type at all the major loci known to influence plumage color. We used pooled sequencing to map the loci affecting feather pigmentation patterns. Our data show that the major locus controlling the segregation of autosomal barring in this cross is located at the distal end of chromosome 11 where the E/MC1R locus is located. We also show that Db is not required for autosomal barring but that it affects the manifestation of this phenotype, and we did not find any evidence for the segregation at the proposed Pg locus.

2 MATERIALS AND METHODS 2.1 Animals

The Fayoumi population used in this study has been kept at an experimental farm over several decades and has shown none or little deviation from the characteristic pigmentation phenotype described for this breed. Feather samples and phenotype data were obtained from a Fayoumi backcross which was produced at the PEAT Poultry Experimental Facility (INRAE, Nouzilly; https://doi.org/10.15454/1.5572326250887292E12). Five silver autosomal barred Fayoumi dams (presumed genotype: Pg/Pg, Db/Db, S/W) were crossed with two Light Brown Leghorn males (wild-type genotype for all feather color loci to the best of our knowledge). From the F1 generation, twelve gold autosomal barred females (presumed genotype Pg/pg+, Db/db+, s+/W) were selected and crossed with another Light Brown Leghorn male to generate 365 backcross progenies in three batches with complete genotype and phenotype data, except that sex recordings were missing for two individuals. The chicks were phenotyped at hatch and at 12 weeks of age as well as genotyped for Db/SOX10 and the E*R(Fay) allele. This protocol received the permit 02410.02 delivered by the French authority for animal experiments, after advice from the INRA Val de Loire ethical committee for animal experimentation.

Feather samples for gene expression were collected after sacrificing the animals, using electronarcosis followed by bleeding. The feather shafts were frozen in liquid nitrogen and stored at −80°C until further processing. DNA samples for genotyping were obtained from blood samples. Blood samples were taken from the wing vein and submerged with EDTA as an anti-clotting agent. Blood samples from other breeds showing various patterning phenotypes were obtained from fancy breeders.

2.2 Whole-genome resequencing of the Fayoumi backcross

For the Fayoumi backcross, between 71 and 102 DNA samples per phenotype category were pooled in equimolar quantities and sequenced to 45× coverage on Illumina HiSeq2500 with 2 × 125 paired-end reads.

Sequencing adapters and low-quality bases were removed using Trim Galore! (with “--stringency 6 -q 15”). The trimmed reads were aligned to the chicken reference genome (GRCg6a/Galgal6) using bwa mem with the “-M” option. PCR duplicates were identified using Picard MarkDuplicates, and variant calling was done using the GATK toolkit (McKenna et al., 2010) v4.1.1.0. First, base quality scores were recalibrated using known variants in dbSNP150. Next, the HaplotypeCaller was used to call variants in gVCF format. All gVCF files were combined and the raw variants were scored and filtered by the VariantRecalibrator module using dbSNP150 as a training set (--truth-sensitivity-filter-level 90.0), followed by a hard filter which excluded SNPs with QD <2.0, FS >60.0, MQ <40.0, MQRankSum < −12.5, or ReadPosRankSum < −8.0 and indels with QD <2.0, FS >200.0 or ReadPosRankSum < −20.0.

Variants used in the linkage analysis were selected based on an earlier version of the analysis, which was performed as described above with the following differences: alignments were done against Galgal4 using bwa aln followed by bwa sampe, the gatk toolkit v3.2.2 was used, an IndelRealigner step was included, all recalibration steps were done using dbSNP140, and “--ts_filter_level 95.0” was used for the VariantRecalibrator. All genomic coordinates from this analysis have been updated to Galgal6 using the liftOver tool.

2.3 Calculation of the fixation index (FST)

FST was estimated for each called variant as FST = (πbetween – πwithin) / πbetween. Variant-level FST values were then averaged across 30 kb sliding windows and visualized as a Manhattan plot after excluding regions with a very low number of variants. The scores for each region were transformed to a Z score, and regions with a score corresponding to p < .05 after Bonferroni correction were reported.

2.4 Genotyping and Sanger sequencing of MC1R and Db/SOX10

The entire Fayoumi pedigree was genotyped for Db, an 8.4 kb deletion upstream of SOX10, by using primers and PCR conditions described elsewhere (Gunnarsson et al., 2011). The MC1R genotype was first evaluated by using an allele discrimination assay on a 7900 HT Fast Real-Time PCR System machine (LifeTechnologies). In brief, 10 ng of DNA was amplified in a reaction using 1× TaqMan Universal PCR mastermix (LifeTechnologies) with 1× Assay mix containing probe and primers and filled up with water to a total reaction volume of 5 μl/sample. The reaction mix was subjected to 40 cycles as follows: 15 s at 92°C and 1 min at 60°C, preceded by 1 × 10 min at 95°C. The data were analyzed using the software sds 2.3. Custom-made probes for the L133Q mutation were obtained from Thermo Fisher.

In order to select genetic markers for fine mapping, we used the resequencing data and extracted SNPs located between position 40,155,407 and 59,843,680 bp on chromosome 1, which occurred at a frequency of 0.4 to 0.6 in the “autosomal barring” pool and <0.2 frequency in the “wild-type” pool. Each SNP had to be covered by at least 20 reads. The same procedure was used for retrieving SNPs from chromosome 11 located between 15,539,481 and 19,932,279. Some high FST SNPs (FST > 0.35) between 6,258,804 bp and 71,048,324 on chromosome 2 were also included. The final composition of SNPs was chosen randomly by keeping the distance between neighboring SNPs in the range 189–322 kb on chromosomes 1 and 11. For chromosome 2, the SNPs were spread out across the chromosome, as there was no distinct haplotype. The genotyping service was provided by Neogen Europe Ltd, Geneseek. The haplotype composition was evaluated using Excel.

2.5 Linkage analysis

Linkage analysis was performed using the cri-map 2.504 software and 76 SNPs obtained through the custom-made genotyping assay as well as genotype data at MC1R and Db/SOX10 (Tables S2 and S3). The initial analyses were carried out using the entire pedigree. Genetic studies have predicted Pg to be located on chromosome 1 (Carefoot, 1999; Moore & Smyth, 1972), but evaluation of phenotype and genotype at Db in our cross suggested that Db is not required for autosomal barring. We therefore eliminated backcross individuals with wild-type plumage and only used those that showed pattern (autosomal barring or unclear) and carried the Fayoumi allele at MC1R as those could reveal a potential-associated locus on chromosome 1 if there was any. We also noticed a strong sex bias within the phenotype group “unclear” (Table 2), which was dominated by females. We therefore performed part of the analysis by only using E*R(Fay)/E*N backcross females.

2.6 Expression analysis and allelic imbalance

RNA extraction from growing feather follicles, primer design, primer testing, gene expression assays as well as allelic imbalance testing was performed as described (Schwochow Thalmann et al., 2017). All primers used for expression analysis are provided in Table S10. Relative gene expression levels of target genes MC1R, NQO1, CDH1, and WWP2 were calculated using the ΔΔCt method and normalized with two housekeeping genes: eukaryotic translation elongation factor (EEF2) and β-actin (Schwochow Thalmann et al., 2017). The following sequence variants were used for the allelic imbalance assays: MC1R (located at chr11:18,841,043), NQO1 (T to C change; chr11:19,037,574), and CDH1 (A to T change; chr11:18,874,575).

2.7 Immunohistochemistry and in situ hybridization

For section immunostaining and in situ hybridization, fixed skin tissue was embedded in paraffin and sectioned at 6–7 µm. After de-paraffination, sections were processed for immunohistochemistry or in situ hybridization. The MITF antibody was from Abcam (ab12039, 1:200 dilution). The peroxidase staining was used after primary antibody treatment as described (Jiang & Chuong, 1992). Non-radioactive in situ hybridization was performed as described (Chuong et al., 1996). Briefly, the sections were treated with proteinase K (10 µg/ml in PBS) for 20 min, re-fixed with 0.2% glutaraldehyde/4% paraformaldehyde, and rinsed with PBT. The sections were then prehybridized in hybridization buffer (containing 50% formamide, 5× sodium citrate/sodium chloride buffer, 1% sodium dodecyl sulfate, 50 µg/ml heparin, 50 µg/ml tRNA) at 65°C for 1 hr. After prehybridization, sections were placed in new prehybridization buffer containing 1–3 µg/ml digoxigenin-labeled riboprobes and hybridized overnight at 65°C. Finally, sections were incubated with alkaline phosphatase-conjugated anti-digoxigenin Fab (Roche, Indianapolis, IN) overnight. Positive signals were detected by incubating the specimens with NBT (nitro-blue tetrazolium)/BCIP (5-bromo-4-chloro-3'-indolyphosphate) substrates (Promega, Madison).

2.8 Whole-genome resequencing of chicken pools representing different within-feather patterning phenotypes

Six samples with pooled DNA from five individuals each of the Brahma, Buttercup, Fayoumi, Hamburg, Plymouth Rock, and Sebright breeds were sequenced to about 10× coverage on Illumina HiSeq 2,500 with 2 × 100 bp paired-end reads. Genome alignment and variant calling were done against Galgal6 using bwa mem. Variants were called with HaplotypeCaller from the gatk toolkit (v3.8) followed by filtering using the VariantRecalibrator and hard filters as described above.

Pooled heterozygosity was calculated as described (Rubin et al., 2010). For each observed SNP, the read count of the major allele, nMAJ, and the minor allele, nMIN, was calculated across all breeds. A pooled heterozygosity, Hp, was calculated as Hp = 2 * ∑(nMAJ) * ∑(nMIN ) / (∑(nMAJ) + ∑(nMIN ))2 across all variants within a 30 kb window and then transformed to a Z score according to ZHp = (Hp − μ) / σ. Negative Z-scores correspond to regions with less than expected heterozygosity, indicating that they may be identical by descent (IDB). Windows with Z-scores less than −4.70 were considered significant, which correspond to a p-value less than .05 after Bonferroni correction. The resulting Z score was visualized in Manhattan plot-style, omitting windows with too few variants.

For the visualization of all haplotypes near MC1R (Figure 5), we extracted all SNPs (including low-confidence ones) ±5kb (n = 105) or ±50kb (n = 1,424) of the MC1R gene.

3 RESULTS 3.1 Segregation of plumage color in a Fayoumi backcross population

We crossed five Fayoumi females (presumed genotype Pg/Pg, Db/Db) displaying the typical black and white autosomal barring pattern with two Light Brown Leghorn males showing no patterning (wild-type at the Pg and Db loci) (Figure 1a). Twelve F1 females (heterozygous carriers of Pg and Db) showing the typical autosomal barring pattern of the Fayoumi breed were backcrossed with another Light Brown Leghorn male generating a total of 365 backcross progeny with complete phenotype information that were used for the genetic analysis. Of these 365, 102 (43 males and 59 females) chickens exhibited the typical autosomal barring pattern (AB), whereas 203 (99 males, 102 females, and 2 unrecorded sex) were classified as wild type (WT) for patterning (Table 1, Figure S2). A total of 60 birds were more difficult to categorize, 28 (9 males and 19 females) were neither plain nor did they have the typical autosomal barring pattern and were classified as “unclear” (Figure 2a). Furthermore, for 21 males initially phenotyped as AB at hatch and 7, which exhibited an irregular pattern, the pattern had completely disappeared at 12 weeks of age (Table S1, Figure S2), and instead, they developed a reddish taint in the belly region and were termed “red belly” (RB). This group also includes four males, which were scored as wild-type at hatch but developed the red belly phenotype by 12 weeks. These 28 + 4 male progenies classified as red belly were excluded from most analyses.

TABLE 1. Segregation of four different phenotypes in the Fayoumi/Light Brown Leghorn backcross Sex Phenotype at 12 weeks of age Total AB WT unclear RB Males 43 99 9 32 183 Females 59 102 19 0 180 Unknowna 0 2 0 0 2 Total 102 203 28 32 365 Abbreviations: AB, autosomal barring; RB, red belly; unclear, irregular pattern; WT, wild-type.

A 1 Mb region on chromosome 11 is associated with autosomal barring. (a) Segregation of pigmentation phenotypes among backcross progeny. Photograph credit: David Gourichon, INRA. (b) Genome-wide screen for genetic differentiation between three different pools (“autosomal barred,” “unclear” and “wild-type”) of Fayoumi/Light Brown Leghorn backcross offspring using FST values in 30 kb sliding windows. (c) SNP scoring within 20 Mb and 3.5 Mb regions on chromosomes 1 and 11, respectively, in parents (P) and offspring progeny (BC) showing autosomal barring (AB), “unclear” patterning or wild-type plumage. The positions of MC1R and SOX10 are indicated. Please note that a higher density of SNPs was used toward the end of chromosome 11 visually suggesting that the 1 Mb interval is taking up half of the 3.5 Mb region. The borders of the 1 Mb non-recombining interval are based on a single recombinant. Yellow color indicates homozygosity for an allele inherited from the wild-type Light Brown Leghorn parental (WT), whereas red indicates homozygosity for a Fayoumi-derived allele. Orange color indicates heterozygosity. The panel for chromosome 11 represents 18 SNPs with the first situated at 15.8 Mb and the last one at 19.7 Mb. Forty-six SNPs were placed within the 20 Mb region on chromosome 1 with the first SNP located at 40 Mb and the last at 59.8 Mb. White fields indicate missing genotypes.

According to the assumed inheritance pattern of autosomal barring, we expected that 25%–50% of the backcross progeny should be heterozygous carriers of Pg and Db and show this phenotype, the reason for this range in expected frequency is because Pg has been reported to be linked to Db but no precise estimate of the recombination rate is available. We observed 102 AB versus 201 WT, 28 birds with unclear phenotype classification, and 32 with the red belly phenotype (Table 1). By genotyping, we confirmed homozygosity for the causal mutation at Db/SOX10 (Gunnarsson et al., 2011) in all Fayoumi founder females (Db/Db) and all F1 females were heterozygous Db/N as expected. Surprisingly, not all of the autosomal barred backcross progeny were genotyped as carriers of Db based on the presence of the deletion upstream of SOX10 (66 heterozygous carriers of Db and 36 non-carriers) whereas the “unclear” group was almost entirely (27 out of 32 chickens) wild-type for the SOX10 deletion (Table 2). The WT group contained an almost equal distribution with 100 heterozygous carriers of Db and 103 non-carriers. Thus, the segregation and genotyping data question the critical role of the SOX10 deletion for autosomal barring as well as the linkage between Pg and Db.

TABLE 2. Phenotype—genotype associations among backcross progenies of the Fayoumi/Light Brown Leghorn cross Genotype Phenotype AB WT unclear RB Total SOX10 Db*N/N 36 103 25 27 191 Db*Db/N 66 100 3 5 174 MC1R E*N/N 0 203 0 3 206 E*R(Fay)/N 102 0 28 29 159 Total 102 203 28 32 365 Abbreviations: AB, autosomal barring; RB, red belly; unclear, irregular pattern; WT, wild-type. 3.2 A major locus underlying autosomal barring maps to a region on chromosome 11

Based on the phenotype data and the genotype data at the Db/SOX10 locus described above, we set up four pools for whole-genome resequencing using Illumina HiSeq technology in an attempt to map the Pg locus. We used two different wild-type pools but with different genotypes at Db/SOX10 (heterozygous carrier of Db, n = 99 or wild-type, n = 102) whereas the two other pools either contained clearly autosomal barred chickens (heterozygous carrier of Db, n = 71) or chickens with a less clear barring phenotype including all “unclear” progenies as well as some with a less pronounced autosomal barring pattern (wild-type at Db, n = 72). The pools were constructed taking into account the genotype at the Db/SOX10 locus in order to maximize the chance to detect other loci affecting the patterning phenotype. All four pools were compared to each other, and regions of high differentiation between pools were determined using the Fixation index (FST) in 30 kb sliding windows (Figure 2b; Figure S3). Depending on which pools were compared, two regions with high FST, one on chromosome 1 and one on chromosome 11, were detected. The signal on chromosome 1 encompassing the Db/SOX10 locus was expected since we constructed the pools on the basis of the genotypes at this locus. The signal on chromosome 11 covered the region 18.399–19.527 Mb (near the chromosome end) with the highest signal around 18.702–19.185 Mb, which included MC1R located at 18,840,646–18,841,590 bp. The signal on chromosome 11 was detected both when the “autosomal barring” pool and the “unclear” pool were compared with either one or both of the wild-type pools (Figure 2b).

3.3 Fine mapping revealed a region of low recombination on chromosome 11

We used our resequencing data to identify a total of 100 highly informative SNPs. The SNPs were selected from the candidate regions on chromosome 1 (n = 63) and chromosome 11 (n = 21), and in addition a few relatively high FST SNPs from chromosome 2 (n = 16) (Table S2). The entire pedigree was genotyped for these markers to validate the observed associations and improve the map resolution on chromosomes 1 and 11 (Table S3). A total of 46 SNPs on chromosome 1 were informative and revealed no haplotype that was shared among all individuals of either phenotype category (Figure 2c). However, a large proportion of the backcross progeny showing autosomal barring carried the Db haplotype whereas only three that were classified as unclear carried this haplotype revealing a clear association to a patterning phenotype. For chromosome 11, 20 SNPs were highly informative and revealed a region of 1 Mb without any recombination starting from about 18.8 Mb to the end of chromosome 11 (Figure 2c, Table S3). All birds phenotypically classified as autosomal barred or unclear were heterozygous for the MC1R haplotype (E*R(Fay)) inherited from the Fayoumi founder females. Linkage analysis revealed highly significant LOD scores for loci located at the distal end of chromosome 11 (LOD scores = 28.4–69.1), whereas only weak associations between autosomal barring and some markers on chromosome 1 were observed, irrespective of which subset of backcross progeny was used (Tables S4–S6).

The 1 Mb interval on chromosome 11 showing a complete association with patterning contains 29 genes including MC1R (Table S7). We used our resequencing data to extract SNPs that differed between autosomal barred and wild-type backcross progeny within the non-recombining region. The L133Q mutation was not extracted, most likely because it did not meet the coverage requirements of at least 20 reads/SNP. The obtained variant SNPs were analyzed using the online tool UCSC Variant Annotation Integrator (VAI), which resulted in a total of 6,627 variants. As expected, the great majority of those variants were located in non-coding regions such as intergenic regions (2,342), introns (2,298) as well as down- or upstream of genes (967 and 852, respectively; Figure S4), 19 variants were detected as potential splice variants and 149 as exonic variants. Of the 149 exonic variants, 107 were synonymous while 41 were non-synonymous (Figure S4). Ten of the genes in the interval showing no recombination harbored non-synonymous changes (Table S8). To further predict the possible effect of the missense mutations on protein function, we used PROVEAN (http://provean.jcvi.org/index.php), an online tool, which generates a score indicating how likely it is that a missense mutation affects protein function (Choi & Chan, 2015; Choi et al., 2012). Except for the known L133Q mutation in the MC1R E*R(Fay) allele, none of the other missense mutations associated with the Fayoumi haplotype from this region was predicted as having a deleterious effect on protein function (Table S8).

Thus, the screen for missense mutations within the 1 Mb region on chromosome 11 revealed only one missense mutation, the one in MC1R, that was predicted to affect protein function. Since MC1R is an obvious candidate gene for a pigmentation phenotype, we genotyped the entire pedigree for this missense (L133Q) mutation (Table 2). All chicken displaying autosomal barring (AB) or some kind of pattern (“unclear”) carried the variant allele at this position, whereas all plain (WT) chicken had the wild-type allele at MC1R. Even the 28 males that initially were phenotyped as AB but did not show any patterning at 12 weeks of age and were grouped into the red belly group carried the variant allele (Table S1). The segregation data at the MC1R locus among the backcross progeny deviated significantly from the expected 1:1 ratio (χ2 = 6.05, df = 1, p < .05). This is not specific for MC1R but applies to all markers in the 1 Mb region that did not recombine in this pedigree data.

To summarize, our genotyping data revealed a perfect association between a patterning phenotype and the variant E*R(Fay) allele at MC1R. It further revealed an incomplete association between the Db allele at SOX10 and autosomal barring, which suggests that Db is not required to form the pattern but makes it more pronounced and easier to phenotype. There may be a third locus segregating in our cross, which is responsible for the red belly phenotype only present in males. Evaluating the genetic basis for this phenotype will be a subject for future studies.

3.4 MC1R, NQO1, and CDH1 are over-expressed in barred feathers but do not show allelic imbalance

We decided to explore possible regulatory changes in gene expression for genes located in the 1 Mb interval on chromosome 11 associated with patterning and with a putative role in melanocyte biology (Figure 3a; Table S7). In addition to the obvious positional candidate MC1R, we considered the NAD(P)H Quinone Dehydrogenase 1 (NQO1) gene that codes for a cytoplasmatic 2-electron reductase, which is reducing quinones to hydroquinones and has been shown to affect the regulation of tyrosinase, thereby enhancing melanogenesis (Choi et al., 2010). We also examined the expression of Cadherin 1 (CDH1), a calcium-dependent cell-cell adhesion protein, which among others, anchors melanocytes to surrounding keratinocytes (Vasioukhin et al., 2000) and is implicated in diseases such as vitiligo (Tarle et al., 2015). The last gene we evaluated is coding for an E3 ubiquitin-protein ligase (ww domain-containing protein 1, WWP1) involved in the tanning response following UV exposure in melanocytes (Cao et al., 2013). We measured the expression of these genes in growing chicken feathers from 10 to 12 non-barred and 10 autosomal barred birds. Three out of the four genes showed a statistically significant up-regulation of expression in barred feathers compared with non-barred feathers (Figure 3b): MC1R (fold change 2.3 ± 0.2, p = .0007, Student's t test), NQO1 (fold change 3.2 ± 0.7, p = .004) and CDH1 (fold change 1.7 ± 0.2, p = .008). WWP1 showed no statistically significant differential expression (p = .75, Student's t test).

Expression analysis of candidate genes on chromosome 11. (a) Genes in the interval lacking recombination in the Fayoumi/Light Brown Leghorn cross. Positional candidate genes are highlighted with red text. (b) Relative gene expression levels of candidate genes on chromosome 11 in autosomal barred and non-barred feathers. EEF2 and βACTIN were used as housekeeping genes. Significant differences are indicated by stars (Student's t test; *p < .05, **p < .01, ***p < .001). (c) Relative proportion of the Fayoumi (F) allele at MC1R, NQO1, and CDH1 in either cDNA samples from F/- feathers (light gray) or genomic DNA from different genotypes (dark gray; F/-, -/- and F/F)

We next reasoned that if the elevated expression levels are the result of cis-regulatory changes, we expect to observe allelic imbalance in gene expression at one or more of the three genes (Figure 3c). We used SNPs in the transcripts to assess the relative expression of alleles and genomic DNA as control in which we expect a perfect 50:50 ratio. We used the T to A polymorphism causing the MC1R/E*R(Fay) allele (L133Q; chr11: 18,841,043 bp), and SNPs in NQO1 (T to C change; chr11: 19,037,574 bp) and CDH1 (A to T change; chr11:18,874,575 bp) to assess allelic expression in 15 heterozygous chickens utilizing cDNA from growing feathers and the pyrosequencing technology. The expression of the Fayoumi allele associated with each investigated gene was 43.1 ± 0.75% for MC1R, 47.7 ± 0.43% for NQO1, and 59.4 ± 1.11% for CDH1. These cDNA data and the genomic DNA control differed significantly from each other (MC1R: Student's t test p = .0001; NQO1: Student's t test p = .008; CDH1: Student's t test p = 1 × 10−6), but these minor differences do not support a typical allelic imbalance expected in the presence of cis-regulatory effects, and we believe that they are unlikely to reflect actual expression differences with biological significance.

To summarize, MC1R, NQO1, and CDH1 located within the non-recombining region show up-regulated expression in growing autosomal barred feathers but this up-regulation does not appear to be mediated by a cis-regulatory effect since none of the three genes shows clear allelic imbalance in favor of the Fayoumi allele. A possible explanation for this up-regulated expression is that the Fayoumi allele underlying autosomal barring results in a higher proliferation of melanocytes.

3.5 The expression of MITF, ASIP, and KIT in Fayoumi feather follicles

Pigment bars can form due to the presence/absence of melanocytes or by differences in melanin production (eumelanin, pheomelanin, or no melanin) of the melanocytes that are present (Andersson et al., 2020; Lin et al., 2013). Autosomal barring is composed of alternating black and non-black bars (Figure 1). The non-black bars may be yellowish due to the presence of pheomelanin (Figure 1b) or more or less non-pigmented, like in Fayoumi chicken (Figure 1a) due to the presence of the Silver allele at SLC25A2 that inhibits expression of pheomelanin (Gunnarsson et al., 2007). To characterize the mechanism underlying autosomal barring in Fayoumi chicken, we performed immunostaining of MITF (microphthalmia-associated transcription factor), a marker of melanocyte progenitor cells. MITF-positive cells are present in the proximal follicle where melanocyte stem cells are present, as shown by positive staining for KIT (Figure 4a,c). Toward the distal feather, barb branches start to form and pigment bars emerge. In both pigmented and non-pigmented regions, we observe the presence of MITF-positive melanocyte progenitors. This is particularly clear in the border region (Figure 4a,b). We explored what may repress the activity of these melanocyte progenitors in the non-black region. In Silver Laced Wyandotte chicken, ASIP was found to be present in the peripheral pulp facing the white region and suppress eumelanogenesis of MITF-positive progenitor cells (Inaba & Chuong, 2020; Lin et al., 2013). Using in situ hybridization, we detected expression of ASIP in the non-black regions. This is particularly clear in a cross section (Figure 4b). Thus, the non-black bars in autosomal barring of Fayoumi chicken are not caused by the absence of melanocytes, as seen in sex-linked barring (Lin et al., 2013; Schwochow Thalmann et al., 2017). Our results suggest that eumelanogenesis is suppressed due to the expression of ASIP in the non-black regions.

The expression of MITF, ASIP, and KIT in the Fayoumi feather follicle. (a) MITF immunostaining. Longitudinal feather sections with enlargement shown in the right column. a–c, black and white barred region. d, proximal follicle in collar bulge region. MITF-positive cells are seen in the basal layer of the feather filament epidermis in the stem cells region, in both black and white barb ridges. MITF-positive cells (red) are highlighted by box a, b, c, and d. (b) Longitudinal feather sections with ASIP in situ hybridization. ASIP is expressed in feather peripheral pulp facing non-pigmented bar regions (a’ and b’; b’ is a cross section). (c) Longitudinal feather sections with KIT in situ hybridization showing KIT-positive melanocyte stem cells at the follicle base near collar bulge region box a’’ and b’’

3.6 Patterning phenotypes are strongly associated with certain MC1R (E) alleles

A number of chicken breeds are assumed to be fixed for the Pg allele based on the interpretation of the genetic basis for various plumage color variants (Figure 1b; Figure S1). We therefore performed whole-genome pooled sequencing of birds from six chicken breeds (n = 5 each), exhibiting a variety of pigmentation patterns: Partridge Plymouth Rock—penciling, Buttercup—autosomal barring, Silver Sebright—single lacing, Brahma—penciling, Hamburg—spangling, and Fayoumi—autosomal barring (Figure 1; Figure S1). This was carried out with two purposes (i) to explore whether all breeds predicted to be fixed for the Pg allele share any region of the genome that is identical by descent (IBD) and (ii) if these breeds share an IBD region within the 1 Mb region on chromosome 11 defined in our Fayoumi backcross population. Thus, the reason for using only five individuals per pool was that our aim was not to estimate allele frequencies, but to identify sequence variants that are fixed in these breeds since all individuals within the breeds show a patterning phenotype. The analysis of the pooled data revealed no striking IBD region shared by all populations neither within the chromosome 11 interval nor in the entire genome (Figure S5). However, the most striking finding was that all populations carried only two alleles Birchen (E*R or E*R(Fay)) or brown (E*B) at MC1R, while not sharing any common haplotype around MC1R (Figure 5).

A striking finding when comparing the 10 kb region harboring MC1R is the lack of linkage disequilibrium between sequence variants in this region (Figure 5), a total contrast to the lack of recombination over a 1 Mb region reported in this study based on pedigree analysis. This analysis also demonstrates that the Birchen (E*R) allele previously defined based on its phenotypic effect on plumage color is genetically heterogenous both as regards the MC1R coding sequence as well as the flanking sequences that may harbor regulatory variants affecting MC1R expression. For instance, the Brahma breed is considered fixed for the E*R allele but the Brahma birds included in this study are apparently segregating for at least two alleles both carry