# A novel locus on chromosome 1 underlies the evolution of a melanic plumage polymorphism in a wild songbird

Yann X. C. Bourgeois, Boris Delahaie, Mathieu Gautier, Emeline Lhuillier, Pierre-Jean G. Malé, Joris A. M. Bertrand, Josselin Cornuault, Kazumasa Wakamatsu, Olivier Bouchez, Claire Mould, et al.

► To cite this version:

Yann X. C. Bourgeois, Boris Delahaie, Mathieu Gautier, Emeline Lhuillier, Pierre-Jean G. Malé, et al. A novel locus on chromosome 1 underlies the evolution of a melanic plumage polymorphism in a wild songbird. Royal Society Open Science, 2017, 4 (2), Non paginé. 10.1098/rsos.160805. hal-01603958

HAL Id: hal-01603958

https://hal.science/hal-01603958

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ROYAL SOCIETY
OPEN SCIENCE

rsos.royalsocietypublishing.org

Research

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

**Cite this article:** Bourgeois YXC *et al*. 2017 A novel locus on chromosome 1 underlies the evolution of a melanic plumage polymorphism in a wild songbird. *R. Soc. open sci.* **4**: 160805.
http://dx.doi.org/10.1098/rsos.160805

Received: 11 October 2016

Accepted: 12 January 2017

**Subject Category:**

Biology (whole organism)

**Subject Areas:**

ecology/evolution/genomics

**Keywords:**

melanism, chromosome 1, polymorphism, selective sweep, *Zosterops*

**Authors for correspondence:**

Yann X. C. Bourgeois

e-mail: yann.x.c.bourgeois@gmail.com

Borja Milá

e-mail: b.mila@csic.es

Christophe Thébaud

e-mail: christophe.thebaud@univ-tlse3.fr

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.fishshare.c.3679084.

THE ROYAL SOCIETY
PUBLISHING

A novel locus on
chromosome 1 underlies
the evolution of a melanic
plumage polymorphism
in a wild songbird

Yann X. C. Bourgeois<fnref n="1" />,†, Boris Delahaie<fnref n="1" />,
Mathieu Gautier<fnref n="2" />, Emeline Lhuillier<fnref n="3" />,<fnref n="4" />,
Pierre-Jean G. Malé<fnref n="1" />, Joris A. M. Bertrand<fnref n="1" />,
Josselin Cornuault<fnref n="1" />, Kazumasa Wakamatsu<fnref n="5" />,
Olivier Bouchez<fnref n="3" />,<fnref n="6" />, Claire Mould<fnref n="1" />, Jade Bruxaux<fnref n="1" />,
Hélène Holota<fnref n="1" />, Borja Milá<fnref n="7" />,‡ and
Christophe Thébaud<fnref n="1" />,‡

<figure><img src="image_4.png" /></figure>

YXCB, 0000-0002-1809-387X; EL, 0000-0003-2629-7773

Understanding the mechanisms responsible for phenotypic diversification within and among species ultimately rests with linking naturally occurring mutations to functionally and ecologically significant traits. Colour polymorphisms are of great interest in this context because discrete colour patterns within a population are often controlled by just a few genes in a common environment. We investigated how and why phenotypic diversity arose and persists in

<footnote>†Present address: New York University Abu Dhabi, PO Box 129188, United Arab Emirates.</footnote>

<footnote>‡These authors contributed equally to this study.</footnote>

<footnote n="1">¹Laboratoire Évolution et Diversité Biologique, UMR5174 CNRS, Université Paul Sabatier – ENFA, 31062 Toulouse Cedex 9, France</footnote>

<footnote n="2">²INRA, UMR 1062 CBGP (INRA, IRD, Cirad, Montpellier SupAgro), Campus de Baillarguet, 34988 Montferrier-sur-Lez, France</footnote>

<footnote n="3">³INRA, GeT-PlaGe, Genotoul, 24 chemin de Borde Rouge, Auzeville, CS 52627, 31326 Castanet-Tolosan, France</footnote>

<footnote n="4">⁴INRA, UAR1209, 24 chemin de Borde Rouge, Auzeville, CS 52627, 31326 Castanet-Tolosan, France</footnote>

<footnote n="5">⁵Department of Chemistry, Fujita Health University, School of Health Sciences, Toyoake Aichi 470-1192, Japan</footnote>

<footnote n="6">⁶GenPhySE, Université de Toulouse, INRA, INPT, INP-ENVT, 24 chemin de Borde Rouge, Auzeville, CS 52627, 31326 Castanet-Tolosan, France</footnote>

<footnote n="7">⁷National Museum of Natural Sciences, Spanish National Research Council (CSIC), 28006 Madrid, Spain</footnote>
the *Zosterops borbonicus* white-eye of Reunion (Mascarene archipelago), a colour polymorphic songbird in which all highland populations contain individuals belonging to either a brown or a grey plumage morph. Using extensive phenotypic and genomic data, we demonstrate that this melanin-based colour polymorphism is controlled by a single locus on chromosome 1 with two large-effect alleles, which was not previously described as affecting hair or feather colour. Differences between colour morphs appear to rely upon complex *cis*-regulatory variation that either prevents the synthesis of pheomelanin in grey feathers, or increases its production in brown ones. We used coalescent analyses to show that, from a 'brown' ancestral population, the dominant 'grey' allele spread quickly once it arose from a new mutation. Since colour morphs are always found in mixture, this implies that the selected allele does not go to fixation, but instead reaches an intermediate frequency, as would be expected under balancing selection.

# 1. Introduction

Identifying the specific genes underlying inter-individual phenotypic variation and reconstructing their evolutionary history is a key issue to link naturally occurring mutations to ecologically significant traits, and help resolve questions relative to the origin and maintenance of genetic and phenotypic variability in natural populations [1]. In this context, colour polymorphisms have played an essential role by enhancing our understanding of how selection and demography can impact phenotypes [2–5], and by identifying proximate causes of phenotypic variation [6,7]. Melanins, in contrast to other pigments such as carotenoids, allow the production of colour traits that in most cases are independent of the environment and the individuals’ phenotypic condition [8,9] (but see also [10,11]) and show heritable segregation among colour morphs [12], a feature that appears to be critical for explaining the origin and maintenance of distinct colour morphs in natural populations [13]. Moreover, the observed covariation between melanin-based pigmentation and life-history or social strategies [14,15] suggests a major role for intergenic interactions and pleiotropic effects in the evolution of colour polymorphisms [16].

A great deal of research over the last decade has been devoted to elucidating the molecular basis of melanin-based colour polymorphism in natural populations. Early studies have assessed the role of the melanocortin-1-receptor (*MC1R*) gene, currently one of the most studied ‘colour genes’. This gene displays relatively low pleiotropic effect (but see [17]) and has been shown to explain colour variation in a broad variety of independent lineages, such as lizards [18], mice [19], humans [20], mammoths [21] and birds [22–25]. Although this suggested at first a relatively simple and general genetic mechanism underlying melanin-based polymorphisms, there is mounting evidence suggesting that mechanisms of hair/feather colour evolution may involve many other genes, and this seems particularly likely in species or populations in which the polymorphism arises from complex patterns of eumelanin/phaeomelanin deposition [26,27]. Studies in model species such as laboratory mice have shown that several genes other than *MC1R* could play a role in explaining differences in the patterning of melanin pigments across the body [28], and alternative candidate loci for melanic colour polymorphisms have indeed been identified in wild non-model species [29–32].

The recent development of next-generation sequencing has allowed research to transcend candidate-gene approaches in wild non-model species, leading to the identification of genomic regions associated with fitness-related traits, including colour traits [33,34]. Here, we took advantage of this technology and used a population genomic approach to investigate the genetic architecture underpinning a melanic plumage polymorphism in natural populations of the Reunion grey white-eye, *Zosterops borbonicus*, a textbook example of intraspecific variation in plumage coloration in birds [35]. This species displays four geographically structured plumage forms distributed across the small oceanic island of Reunion (2512 km²) [36]. One of these forms, restricted to the highlands of Reunion, comprises two distinct and sympatric colour morphs, with birds showing predominantly grey or brown plumage, respectively (figure 1). Melanic plumage polymorphism is widely maintained across the range of this highland form, with both grey and brown birds always present at any given locality, although populations vary in morph frequencies [30]. Field observations suggesting that mating between grey individuals could produce both grey and brown offspring, and that mating between brown individuals always produced brown offspring, led to the proposition that differences between grey and brown morphs could be due to genetic changes at a few loci of major effect [30].

Previously, we used a candidate-gene approach to search for associations between plumage colour and genetic variants in the Reunion grey white-eye, and did not detect any for a series of genes
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Phenotypic variation in highland populations of the Reunion grey white-eye and the population sampling scheme. (a) Photographs of brown and grey morphs (i), microscopic view of grey and brown feathers (ii) and melanin content per morph (iii) for 12 brown and 12 grey individuals. (b) PCA on spectrophotometry measures for feathers sampled from the back of brown and grey birds. (c) Morph frequencies across different localities. 'P' stands for localities included in pedigree analyses. Localities used for GBS and RAD-seq analyses are labelled. The size of the pie charts is proportional to sample size (ranging from 11 to 111 birds).</figcaption></figure>

<table><caption>**Table 1.** Summary of datasets used in this study.</caption><thead><tr><td>method used</td><td>goal</td><td>notes</td><td>sample size</td><td>references</td></tr></thead><tbody><tr><td>microsatellites</td><td>pedigree analysis to characterize inheritance patterns of plumage colour</td><td>obtained from two populations monitored since 2008 (Pas de Bellecombe and Bois Ozoux) and one population sampled in 2007 and 2012 (figure 1)</td><td>261</td><td>[39]</td></tr><tr><td>pooled RAD-sequencing</td><td>mapping loci associated with plumage colour</td><td>high density of markers (more than 600 000)</td><td>137 (6 pools)</td><td>[40]</td></tr><tr><td>GBS</td><td>validating associations obtained from pooled data with individual genotyping</td><td>lower density of markers than RAD-seq (25 000)</td><td>42</td><td>this study</td></tr><tr><td>whole-genome resequencing</td><td>annotating coding/non-coding mutations; characterizing selection in the genomic region associated with plumage colour</td><td>near-exhaustive coverage of the region associated with plumage colour</td><td>12 (6 parents and 6 offspring)</td><td>this study</td></tr></tbody></table>
that play a critical role in melanin-based pigmentation patterns in other species, including *MC1R*, *Agouti*, *Corin*, *Pro-opiomelanocortin* and *Tyrosinase-related protein 1* [37,38]. In this study, we first confirmed that segregating phenotypes conform to a Mendelian model using genetic marker-based pedigree reconstruction (table 1). Then, using restriction site-associated DNA (RAD)-sequencing [41] on pools of individuals from three geographically separate populations [40], we performed a genome-wide scan relative to colour morphs while accounting for population structure. This approach provides a cost-effective way to detect associations by allowing large numbers of individuals to be included into a sample and by providing a sufficient density of markers for mapping (see table 1 for more details on the approaches used in this study). This led us to discover a single narrow region on chromosome 1 showing a strong association with melanin-based plumage colour phenotypes. We confirmed the existence of this association by using individual genotyping-by-sequencing (GBS) [42] and whole-genome sequencing of individuals with known pedigrees, and assessed the functional consequences of the molecular variants underlying plumage colour variation. Finally, we addressed the origin of the polymorphism by estimating the relative age of the genetic variants and tested the role of selection in explaining how the polymorphism could have spread and been widely maintained across populations.

## 2. Material and methods

### 2.1. Field sampling

Birds were captured in the field using mist nets between 2007 and 2012 on Reunion (55°39' E; 21°00' S), were weighed and marked with a uniquely numbered aluminium ring, and approximately 10 µl of blood was collected from each bird. Blood was conserved in Queen's lysis buffer [43] and stored at -20°C for long-term preservation. Individuals were sexed by PCR [44] in order to infer the number of distinct Z chromosomes included in each pool. On each bird, 10 feathers from four different body parts (head, back, flank and belly) were collected for spectrophotometric analyses. We also measured five morphological traits with a dial calliper (to the nearest 0.1 mm): tail length (from the uropygial gland to the tip of the longest rectrix), tarsus length (from the intertarsal joint to the most distal undivided scute on the tarsometatarsus), bill length (from the anterior end of the nares to the tip of the upper mandible), bill width and depth (both measured at the anterior end of the nares).

### 2.2. Measurement of the reflectance spectra of plumage patches

Ninety-six individuals from the 137 used for pooled RAD-sequencing were used for characterizing phenotypes by spectrophotometry analyses. We used reflectance spectrophotometry to characterize variation in plumage colour, because it provides an objective quantification of colour [45,46]. We summarized colour variation by conducting a principal components analysis (PCA) on a set of five variables summarizing the spectrum for back feathers that were then used for estimating melanin content. We determined whether birds were actually able to discriminate morphs by using Vorobeyv & Osorio's model of colour perception [47]. This model allows chromatic contrasts between two colours to be compared with a threshold value of colour discrimination by birds (for further details of how this was perfomed, see [48]).

### 2.3. Melanin content

We used microanalytical methods to quantify the eumelanin and pheomelanin content of the feathers from brown and grey morphs (see [49] for a detailed description of the method). Feather samples (ca 7 mg) from the back of 12 brown and 12 grey individuals were homogenized with a Ten-Broeck homogenizer at a concentration of 10 mg ml⁻¹ and 100 µl aliquots were subjected to Soluene-350 solubilization [50], alkaline hydrogen peroxide oxidation [51] and hydriodic acid hydrolysis [52]. Values are from single determinations. A500/mg and A650/mg values were substracted by background values of 0.019 and 0.001.

### 2.4. Pedigree analysis

Performing controlled matings in wild populations is challenging, and hence to study inheritance patterns, we inferred pedigrees in three natural populations monitored over multiple years at three
localities. Given the reduced dispersal in *Z. borbonicus* [53,54], samples from these populations were more likely to include parents and their offspring. Pedigree reconstruction was performed using the Bayesian parentage assignment algorithm implemented in the R package *MASTERBAYES* v. 2.50 [55] and genotype data from an informative panel of 11 neutral microsatellites [39]. A total of 261 birds of known sex were genotyped and included in this analysis. When known, a birth cohort identifier was specified to minimize the possibility that individuals appear as potential parents of offspring from the same or a previous cohort. Since *MASTERBAYES* allows the simultaneous use of genetic and phenotypic data that may inform on parentage, the reconstruction was performed twice: a first time using genotypic information only and a second time including both genetic and morphometric data (including body mass). The genotyping error rates for each locus, the number of unsampled sires and the number of unsampled females were all estimated jointly from the pedigree. Markov chains were run for 1.1 million iterations, with a burn-in of 100 000 iterations and a thinning interval of 1000. Parents (or father or mother) were assigned to an offspring with a 95% joint posterior probability threshold.

## 2.5. RAD-sequencing using pooled DNA samples

To identify loci associated with colour, we used a paired-end RAD-sequencing protocol, using a dataset described previously in which six pools of 18–25 individuals representing the two colour morphs in each of three separate populations (named ‘Bois Ozoux’, ‘Pas de Bellecombe’ and ‘Tévelave’) were sequenced [40]. We took advantage of the recent sequencing of the *Z. lateralis* genome [56] to map reads back onto this reference with BWA MEM (v. 0.7.12) [57], instead of creating consensuses directly from data as in [40]. We aligned contigs and scaffolds from *Z. lateralis* on the zebra finch genome (version July 2008, assembly WUGSC v. 3.2.4) using LASTZ [58]. We removed PCR duplicates using samtools [59]. SNPs were called using Popoolation2 (v. 1.201) [60].

## 2.6. Genotyping by sequencing using individual DNA samples

We further validated the results obtained by the pooled RAD-seq approach by performing a genotyping by sequencing analysis [42] on 42 individuals from the same populations that were used to build the pooled DNA samples (seven individuals per morph per population for a total of 14 individuals in each population). Approximately 1 µg of DNA was extracted with a QIAGEN Blood and Tissue kit following the manufacturer’s instructions and sent to the BRC Genomic Diversity Facility at Cornell University [42]. Reads were trimmed with Trimmomatic (v0.33) and mapped on *Z. lateralis* genome using BWA. SNPs were called using freebayes (v0.9.15-1) and filtered with VCFTOOLS (0.1.12b) using the following criteria: (i) a mean sequencing depth between 6 and 20×; (ii) a minimal genotype quality of 20; (iii) a minor allele count of 3, which resulted in around 25 000 unambiguous SNPs.

## 2.7. Genetic structure and association analysis

To validate the lack of genetic structure according to morphs, we performed a PCA [61], using custom scripts in R (v3.2.2) for allele frequencies obtained from pooled data, and the Bioconductor package SeqVarTools for GBS data [62]. Given our sampling strategy, we expected morphs to group together in the PCA and population structure, if any, to correlate with geography. To further rule out the possibility of any neutral structure according to colour morph, we also conducted a locus by locus analysis of molecular variance (AMOVA) in ARLEQUIN v. 3.5 [63] on the 10 115 autosomal SNPs with less than 95% of missing data.

To detect loci displaying a significant association with coloration, we performed an association analysis on the pooled RAD-seq data using the software BAYPASS (v. 2.1) [64]. We computed the empirical Bayesian *p*-value (eBPis) and Bayes factors (BF) expressed in decibans units (dB) to determine the level of association of each SNP with the grey/brown trait. Decibans are a commonly used statistic that describes here the probability that data are produced under a given model. They are often used as a unit for BF. Here, it describes the odds of association to morphs versus the null hypothesis of non-association. BAYPASS was run using default parameters under the core model. Empirical Bayesian *p*-values and BF were then computed for a pseudo-observed dataset containing 1 000 000 SNPs simulated from the actual data. We then compared the observed values to this distribution in order to calibrate the statistics.

For GBS data, we performed a single analysis in PLINK [65] by correcting for population structure with a Cochran–Mantel–Haenszel test (option –mh in PLINK).
## 2.8. Whole-genome sequencing

To validate the hereditary transmission of plumage colour and confirm the validity of pedigrees, 12 individuals from three distinct families were selected for 150 bp paired-end Illumina sequencing. Individuals included 4 brown and 8 grey birds (for reconstructed pedigrees see the electronic supplementary material, table S4). Three types of matings were considered: two grey parents with only grey offspring, two grey parents with both brown and grey offspring, and one parent of each morph with both types of offspring. DNA sequencing libraries were prepared according to Illumina's protocols using the Illumina TruSeq Nano DNA LT Library Prep Kit. Briefly, DNA was fragmented by sonication (Covaris M220) and adaptors were ligated to be sequenced. Eight cycles of PCR were then applied to amplify the libraries. Library quality was assessed using an Agilent Bioanalyzer and libraries were quantified by QPCR using the Kapa Library Quantification Kit. Sequencing was performed in paired-end (2 × 125 pb) on an Illumina HiSeq2500 sequencer at the GeT-PlaGe core facility (INRA, Toulouse). The mean sequencing depth ranged from 6.9× to 11.6× after removing duplicates (9.3× in mean, electronic supplementary material, table S6). After calling with freebayes, SNPs were filtered with VCFTOOLS using the following criteria: (i) a mean sequencing depth between 6 and 20×; (ii) a minimal Phred score of 20 and a genotype quality above 20; (iii) at least eight individuals genotyped.

Putative pedigrees obtained from microsatellite data were confirmed by checking the coefficients of relatedness between individuals in KING (v. 1.4) [66] using 10 000 biallelic SNPs randomly sampled across autosomal scaffolds. An association analysis (including indels) taking into account family structure was performed in LAMP (v. 0.0.12) [67]. We present here results taking into account the putative heredity of the colour trait (brown recessive, grey dominant) as they were similar to those obtained when not taking the transmission mode into account. To assess whether any associated SNP might lead to changes in protein sequence, we performed an annotation using the software SNPdat (v. 1.0.5) [68] and using gene coordinates from [56].

## 2.9. Detecting selective sweeps and characterizing derived and ancestral states

To determine when and on which specific allele selection occurred, we used the coalescent framework implemented in ARGWeaver [69]. ARGWeaver models the coalescent process across non-recombining blocks of sequences and thus provides access to the evolutionary history of DNA sequences. It allows recovering several statistics that describe local genealogies, like coalescence times and local effective population sizes. In the case of a recent and partial selective sweep, selected lineages should display shorter coalescence times than the ancestral ones, i.e. alleles under selection will tend to be younger than neutral alleles.

We ran the analyses using whole-genome sequences from the three parental pairs, as these individuals displayed low levels of relatedness in the whole-genome association analysis and were predicted to display an equal number of brown and grey alleles. To reduce computational burden, we performed the analysis on the scaffold covering the candidate region (scaffold 40) and one on chromosome 2 (scaffold 30) with no SNPs associated with colour. The algorithm was run for 1000 iterations. We then extracted trees and half-time TMRCAs from the output.

# 3. Results

## 3.1. Differences in melanin-based pigmentation patterns between morphs are functionally significant

Brown and grey birds were different in coloration when examined in an avian-appropriate colour space, clustering into two distinct groups along the first component axis in a PCA (figure 1; electronic supplementary material, table S1). All colour variables differed significantly between grey and brown morphs (electronic supplementary material, table S1). In addition, we found that, for all patches, the differences between morphs are greater than the discrimination threshold in the avian visual space (electronic supplementary material, table S1). Taken together, these results highlight the functional significance in relation to colour perception of the differences in pigmentation patterns that distinguish grey and brown birds.

We analysed the melanic content of feathers by determining eumelanin and pheomelanin concentrations (figure 1; electronic supplementary material, table S2). Although levels of eumelanin were
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Population structure and association study. (a) PCA on allele frequencies for pooled RAD-seq (top) and individual GBS (bottom) data. (b) Genome-wide association analysis performed with BAYPASS for 627 795 RAD-seq loci. Horizontal dashed lines correspond to 0.001% (grey) and 0.01 (black) significance thresholds. (c) Density of highly associated SNPs on scaffold 40 from the *Zosterops lateralis* reference genome. Upper stars: SNPs displaying both highly significant BF and empirical Bayesian *p*-values (above the 0.001% threshold). Lower stars: 7 GBS SNPs associated with a *p*-value < 1 × 10⁻⁴. The region including all 100 kb intervals with at least three highly associated SNPs is shown in orange, and the three intervals with the highest density of associated SNPs are shown in red. LASTZ alignment between scaffold 40 and the zebra finch chromosome 1 is also illustrated.</figcaption></figure>

similar in both morphs, levels of pheomelanin were clearly lower in the grey morph when compared with the brown morph (ANOVA, F₁,₂₀ = 132.817, *p* = 2.77 × 10⁻¹⁰). Pheomelanin was mostly concentrated in feather barbs and rachis (figure 1), suggesting that it is produced and deposited in a timely fashion during feather growth [70]. We also found a highly significant relationship between PC1 scores and the pheomelanin marker (Spearman's *ρ* = 0.834, *p*-value = 7.6 × 10⁻⁶), consistent with the fact that most of the variation in reflectance spectra between morphs is related to variation in pheomelanin content.

## 3.2. Lack of population structure and Mendelian inheritance reveal a true genetic polymorphism with single-locus control

Previous field observations suggested a lack of assortative mating relative to plumage colour in white-eye populations, which should lead to a lack of genetic structure according to morph [30]. Pooled RAD-sequencing and GBS data on 137 and 42 individuals, respectively, from three geographically separate populations revealed that genetic structure was consistent with geography but not with colour morphs. PCA on allele frequencies from pooled RAD-seq data highlighted a clear grouping of brown and grey pools by locality (figure 2a). This lack of structure was also found for individual GBS data (figure 2b) and was further confirmed by an AMOVA (electronic supplementary material, table S3) which revealed a significant effect of population structure in shaping neutral genetic diversity (*F*<sub>st</sub> = 0.022, *p*-value < 1 × 10⁻⁵), yet no effect of colour morph (*F*<sub>ct</sub> = -0.011, *p*-value = 1).

To study the pattern of inheritance of the two colour morphs, we generated pedigrees for 36 putative parent-offspring triads based on the genotyping of 11 microsatellite markers in a total of 260 individuals.
<table><caption>**Table 2.** Link between parent and offspring *Zosterops borbonicus* phenotypes. Results were obtained from a pedigree analysis based on 11 microsatellites and morphometric measurements other than colour. Values correspond to pedigree reconstructions congruent between analyses including and excluding morphometric measurements. +1: found in analysis with genetic data only.</caption><thead><tr><th>parental phenotype</th><th>brown offspring</th><th>grey offspring</th></tr></thead><tbody><tr><td>brown × brown</td><td>16</td><td>0</td></tr><tr><td>grey × grey</td><td>1</td><td>8 (+1)</td></tr><tr><td>grey × brown</td><td>6</td><td>5</td></tr></tbody></table>

Matings involving two brown parents systematically led to brown offspring, while grey offspring always had at least one grey parent (table 2), suggesting a single-locus control with grey dominant to brown.

### 3.3. Genome-wide association analysis identifies a single genomic region associated with colour phenotype

To identify the genomic regions responsible for the observed colour polymorphism, we performed a genome-wide association analysis on a total of 627 795 SNPs obtained from the pooled RAD-sequencing data, and mapped it onto the *Z. lateralis* genome (figure 2b). Thirty-three SNPs had an empirical Bayesian *p*-value (eBPis) higher than 1.48 (at the 0.001% threshold). Among these SNPs, 28 were found on scaffold 40 on the reference genome for *Z. lateralis* which covers approximately positions 15–21 Mb on chromosome 1 of the better assembled zebra finch reference genome (figure 2c). Density of strongly associated SNPs (with an eBPis above the 0.01% threshold) was particularly high around 3 Mb of this scaffold, which corresponds approximately to position 18 Mb on the zebra finch chromosome 1. This pattern was further confirmed by individual GBS data, with which seven variants were found to be associated to coloration with an uncorrected *p*-value between 9.1 × 10⁻⁵ and 8.6 × 10⁻⁷ on scaffold 40, while only two other variants were found to be associated with similar levels of significance on scaffold 172, which is directly upstream of scaffold 40 according to the zebra finch reference genome.

In sum, we found a single genomic region displaying a clear signal of association with phenotype. Within this region, we identified seven candidate genes based on their known role or on the role of homologues in melanic pigmentation (table 3). Among those genes, four—*AP1S2*, *GPM6B*, *RAB9A* and *TRAPPC2*—lie within regions where associated SNP density is highest (figure 2).

### 3.4. Functionally important changes lie outside candidate gene coding regions

Having pinpointed the locus responsible for plumage colour variation, we searched for functional polymorphism at non-synonymous sites in the candidate region by obtaining whole-genome sequences from 12 individuals sampled across three distinct families from different populations (electronic supplementary material, table S4). Coefficients of kinship estimated from whole-genome data for these individuals were consistent with the pedigrees estimated from microsatellite data (electronic supplementary material, table S5). The highest density of associated SNPs was found between 17.7 and 18.7 Mb on the zebra finch chromosome 1 (figure 3a). Only one non-synonymous substitution at a level of significance less than 1 × 10⁻³ was found in *ACE2*, an angiotensin (figure 3a). The induced change led to the substitution of methionin to isoleucin, two hydrophobic amino acids which are extremely similar. Therefore, there is no non-synonymous change that would make a suitable candidate for the observed phenotypic variation. This suggests a role for regulatory mutations in determining the observed phenotype.

### 3.5. Evidence for a selective sweep suggests recent selection on the grey allele

To assess whether selection had a role in the emergence of the observed polymorphism and to identify which of the two allele categories (brown or grey) is ancestral, we sampled ancestral recombination graphs for the entire candidate scaffold. Here, we took advantage of the fact that the individuals used in this test displayed six ‘brown’ and six ‘grey’ alleles at the most strongly associated markers. We found that the candidate region displayed a clear signature of a recent selective sweep, with an extremely short half-TMRCA between positions 2 700 000 and 3 100 000 on scaffold 40 (figure 3b). An examination of the local trees for regions containing SNPs associated with phenotype clearly showed that the most recent
<table><caption>**Table 3.** Summary of seven candidate genes for colour variation in *Zosterops borbonicus*. Information based on OMIM (Online Mendelian Inheritance in Man), a database reporting large-scale genotype–phenotype associations in humans and laboratory mice. Positions correspond to coordinates on zebra finch reference genome.</caption><thead><tr><th>candidate gene</th><th>complete name</th><th>position (from first to last exon)</th><th>role in melanocyte</th><th>impact on phenotype</th><th>references</th></tr></thead><tbody><tr><td>RS1</td><td>retinoschisin</td><td>15 288 297–15 291 738</td><td>retinal cells adhesion, cell–cell interaction</td><td>depigmentation in retinal pigment epithelium</td><td>[71]</td></tr><tr><td>AP1S2</td><td>adaptor-related protein complex 1, sigma 2 subunit</td><td>16 724 833–16 754 005</td><td>AP-1 complex is involved in melanosome genesis and is necessary for TYRP1 to reach the melanosome and produce eumelanin</td><td>mutations on another complex, AP-3, lead to the pearl (Ap3b1; coat hypopigmentation) and mocha (Ap3d1; coat colour dilution) phenotypes in mice</td><td>[72–74]</td></tr><tr><td>GPM6B</td><td>glycoprotein M6B</td><td>17 812 158–17 843 054</td><td>membrane protein involved in neuronal tissues. Up-regulated by MITF</td><td>—</td><td>[75]</td></tr><tr><td>TRAPPC2</td><td>trafficking protein particle complex, subunit 2</td><td>17 877 165–17 880 771</td><td>TRAPP complex is involved in vesicle transport and tethering</td><td>in mice, mutations on the subunit TRAPPC6A lead to hypopigmented patches in the coat and retinal epithelium</td><td>[76]</td></tr><tr><td>RAB9A</td><td>Ras-associated protein 9A</td><td>17 886 048–17 886 836</td><td>regulation of vesicular trafficking. Interacts with BLOC-3, involved in Hermansky-Pudlak syndrome</td><td>Rab38: chocolate phenotype in mice. Rab27a: Griscelli syndrome (hypomelanosis and neurological defects)</td><td>[28,77,78]</td></tr><tr><td>APXL</td><td>apical protein of Xenopus-like</td><td>19 986 240–20 035 565</td><td>melanosome biogenesis and transport. Activity requires the G-protein RAB27A</td><td>ocular albinism</td><td>[79]</td></tr><tr><td>OA1</td><td>ocular albinism1</td><td>20 136 404–20 149 393</td><td>melanosome transport and interactions with cytoskeleton</td><td>ocular albinism</td><td>[80,81]</td></tr></tbody></table>

haplotypes were found in grey individuals only (figure 3c). These haplotypes also harboured the allele associated with grey coloration in the association analysis. Median half-TMRCA was estimated at 1733 generations in the region, significantly shorter than the estimate for the remaining regions included in the analysis (median: 54 460 generations; Wilcoxon rank test, $p < 2.2 \times 10^{-16}$). This indicates that grey coloration has been positively selected and that brown and grey phenotypes represent ancestral and derived conditions, respectively.

## 4. Discussion

We confirmed that highland populations of the Reunion grey white-eye present a genetic colour polymorphism that results from the differential deposition of pheomelanin on feathers, leading to predominantly grey or brown plumage, and is controlled by a simple genetic mechanism with a recessive ‘brown’ allele and a dominant ‘grey’ allele. A genome-wide association analysis designed to search for loci involved in the polymorphism led to the identification of a single genomic region on chromosome 1 characterized by SNPs strongly associated with colour phenotype. A coalescent analysis of selection pinpoints the same region, and suggests recent selection on the grey allele. Annotation of these SNPs
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Test for selection using whole-genome sequencing. (a) Association analysis of whole-genome sequencing data, including four brown and eight grey individuals from three different families. Synonymous and non-synonymous SNPs are highlighted in blue and red, respectively. Stars indicate the positions of the seven candidate genes listed in table 3. Density of SNPs associated with a p-value < 1 × 10⁻³ is also provided. (b) Plot of half-TMRCA (in number of generations) for scaffold 40. Stars indicate the position of the three non-recombining blocks examined in (c). (c) Genealogies observed at three distinct points from the colour locus. All three blocks of sequences included SNPs associated with colour with a p-value < 1 × 10⁻³ in the whole-genome association analysis. Branch lengths represent time in generations. Individual phenotype is indicated in brackets. Individuals 1434 (grey morph) and 314 (brown morph) are found to be homozygous at SNPs strongly associated with colour.</figcaption></figure>

using the *Z. lateralis* reference genome [56] did not reveal any non-synonymous mutations or deletions in coding parts of the genomic region associated with variation between colour morphs, indicating that functionally important nucleotide changes in the Reunion grey white-eye must lie outside coding regions. This suggests a role for *cis*-regulatory mutations that either prevent the synthesis of pheomelanin in grey feathers, or increase its production in brown ones. The first mechanism seems more likely in the present case since the grey phenotype seems to represent the derived condition. The fact that eumelanin is found in both brown and grey feathers, while pheomelanin is only present in brown feathers suggests that the underlying genes may be involved in melanogenesis and melanosome transport.

Importantly, the genomic scaffold associated with plumage colour does not include, to the best of our knowledge, genes previously known to be involved in hair or feather coloration in other species, even in model species such as laboratory mice. Three genes found in this region, *RS1*, *APXL* and *OA1*, are known to affect melanogenesis and pigmentation in retinal pigment cells. In addition to these genes, the colour locus includes genes (*GPM6B*, *RAB9A*, *TRAPPC2*, *AP1S2*) known to be involved in intracellular trafficking and in shaping melanosomes in melanocytes (table 3). Among them, *GPM6B* is of particular interest as it is regulated by *MITF*, a key transcription factor involved in melanocyte specification [75]. However, while changes in *GPM6B* expression were found in a comparison of black and hooded crows [30], there is currently no direct evidence for a role of this gene in feather coloration. Our results also confirm previous findings suggesting that variation at genes often studied in natural populations, such as *MC1R*, *Agouti* or *POMC* [32,82,83], is not directly related to colour variation in the Reunion grey white-eye [38] and provide an impressive illustration of the diversity of the mechanisms underlying melanin-based plumage colour evolution in birds and perhaps other vertebrates. Future
work should focus on gene expression differentiation to define the developmental basis of this colour polymorphism and to determine, for example, if the candidate loci found on chromosome 1 are epistatic over other genes classically found associated with melanin-based colour variation.

A large-effect locus can break up into many small- to moderate-effect mutations [32,84]. Thus, although the pattern of inheritance of plumage colour and the existence of a single strongly associated genomic region points towards a relatively simple genetic mechanism, it remains possible that several of the candidate genes we identified are simultaneously involved in shaping the observed colour variation in this study. If mutations arose independently at several genes involved in melanogenesis, their physical proximity would favour linkage between the best combinations of mutations by sieve. Structural rearrangements such as large inversions covering the genomic region associated to colour, if any, would also stabilize associations between several alleles [85]. In that case, the colour locus that we identified in this study may be a further example of the so-called ‘largesse of the genome’ [86,87]—i.e. the predisposition of certain genomic regions to mediate integrated phenotypic shifts—in vertebrates.

Genes like RS1, APXL or OA1 have been studied in the context of ocular diseases or melanomas [88], in a way that suggests important functions in melanocyte development (table 3). Therefore, they probably play a role in ocular development and in colour recognition, which might have an impact on how individuals recognize and favour conspecifics. In addition, GPM6B is known to be involved in neural development, serotonin uptake and bone formation [89,90]. Thus, this gene could influence feather colour, brain development and behaviour. This feature is particularly interesting, since changes in behavioural strategies associated with visual signalling can explain the persistence of morphs over time [91].

Our work reveals a new genomic region not previously associated with melanic coloration in vertebrates, and underscores the importance of genome-based research on non-model species to understand the genetic basis of ecologically significant traits and their role in phenotypic diversification [92]. Our results suggest a strong selective advantage for the dominant ‘grey’ allele once it arose from a new mutation, leading to its fast spread across all highland populations of the Reunion grey white-eye [93]. Similar patterns of selection on *de novo* mutations at colour genes have been recently described in deer mice (*Peromyscus maniculatus*), where multiple independent mutations were selected for cryptic coloration after the colonization of a novel selective environment [31,32]. In white-eyes, both colour morphs are found at an appreciable frequency in all populations, implying that the selected allele does not go to fixation, but instead reaches an intermediate frequency, as would be expected under balancing selection. This may happen as the result of, for example, heterozygote advantage or negative frequency-dependent selection [94]. Distinguishing between these possibilities is notoriously difficult, and future studies should estimate the relative survival and fertility rates of the different ‘colour’ genotypes in order to determine the nature of selection acting on this colour polymorphism and to come closer to an explanation about its persistence.

**Ethics.** Measurements and manipulations were performed by T.D. and B.M. Manipulation of live birds for ringing, measurements and blood sampling was done under a ringing permit (number 14 631) issued by CRBPO Museum National d’Histoire Naturelle (Paris, France), and in compliance with ethical guidelines at Université Paul Sabatier. A special authorization to work on protected species was issued by Direction Régionale de l’Environnement (DIREN Réunion).

**Data accessibility.** All data are deposited on Dryad (phenotypic data, GBS, http://dx.doi.org/10.5061/dryad.1q9g7) [95] and European Nucleotide Archive repositories (raw reads for whole genomes and pooled RAD-seq, accession nos. ERP020509 and ERP002555). Methods and options are further detailed in the electronic supplementary material. Scripts used in this study are freely available on GitHub at https://github.com/YannBourgeois/Scripts_Pool_RADseq.

**Authors’ contributions.** B.M. and C.T. initiated, coordinated and supervised the project, and coordinated the collection of field samples; Y.B., B.M. and C.T. conceived the study and designed the experiments. Molecular data were obtained by Y.B. with contributions from E.L., J.A.M.B., B.D., J.B., C.M., H.H. and O.B.; pigment data were obtained by K.W.; Y.B. and M.G. analysed the data, with contributions from P.J.M., B.D. and J.C.; and Y.B., B.M. and C.T. wrote the paper with comments from other authors. All authors gave their final approval for publication.

**Competing interests.** We declare we have no competing interests.

**Funding.** This work was supported by Fondation pour la Recherche sur la Biodiversité (FRB), Agence Française pour le Développement (AFD), Agence Nationale de la Recherche (ANR-2006-BDIV002), Centre National de la Recherche Scientifique (CNRS) through a PEPS grant, the National Geographic Society and the ‘Laboratoire d’Excellence’ TULIP (ANR-10-LABX-41). The first author was supported by a MESR (Ministère de l’Enseignement Supérieur et de la Recherche) PhD scholarship during this study.

**Acknowledgements.** We thank Alexandre Roulin, Dieter Ebert, Anne Roulin, Marie Manceau and one anonymous reviewer for insightful comments and assistance that greatly improved the manuscript. We want to thank the
Reunion National Park for granting us permission to conduct fieldwork and to collect blood samples. Thomas Duval, Guillaume Gélinaud, Philipp Heeb, Dominique Strasberg, Ben Warren and Juli Broggi assisted with fieldwork. Isabelle Bourgeois provided assistance with drawing figures.

## References

1. Hoekstra HE. 2006 Genetics, development and evolution of adaptive pigmentation in vertebrates. *Heredity* **97**, 222–234. (doi:10.1038/sj.hdy.6800861)

2. Huxley J. 1955 Morphism in Birds. In *Acta XI International Congress of Ornithology* (eds A Portmann, E Sutter), pp. 309–328. Basel, Switzerland: Birkhäuser Verlag.

3. Ford EB. 1975 *Ecological genetics*, 4th edn. London, UK: Chapman and Hall.

4. Wright S. 1943 An analysis of local variability of flower color in Linanthus parryae. *Genetics* **28**, 139–156.

5. Gray SM, McKinnon JS. 2007 Linking color polymorphism maintenance and speciation. *Trends Ecol. Evol.* **22**, 71–79. (doi:10.1016/j.tree.2006.10.005)

6. Lewontin R. 1974 *The genetic basis of evolutionary change*. New York, NY: Columbia University Press.

7. Manceau M, Domingues VS, Linnen CR, Rosenblum EB, Hoekstra HE. 2010 Convergence in pigmentation at multiple levels: mutations, genes and function. *Phil. Trans. R. Soc. B* **365**, 2439–2450. (doi:10.1098/rstb.2010.0104)

8. McGraw KJ, Hill GE. 2000 Differential effects of endoparasitism on the expression of carotenoid- and melanin-based ornamental coloration. *Proc. R. Soc. Lond. B* **267**, 1525–1531. (doi:10.1098/rspb. 2000.1174)

9. McGraw KJ, Mackillop EA, Dale J, Hauber ME. 2002 Different colors reveal different information: how nutritional stress affects the expression of melanin- and structurally based ornamental plumage. *J. Exp. Biol.* **205**, 3747–3755. (doi:10.1007/s00265-010- 1135-5)

10. McGraw KJ. 2008 An update on the honesty of melanin-based color signals in birds. *Pigment Cell Melanoma Res.* **21**, 133–138. (doi:10.1111/j.1755- 148X.2008.00454.x)

11. Roulin A. 2016 Condition-dependence, pleiotropy and the handicap principle of sexual selection in melanin-based colouration. *Biol. Rev. Evl.* **91**, 328–348. (doi:10.1111/brv.12171)

12. Roulin A, Ducrest A-L. 2013 Genetics of colouration in birds. *Semin. Cell Dev. Biol.* **24**, 594–608. (doi:10.1016/j.semcdb.2013.05.005)

13. Roulin A. 2004 The evolution, maintenance and adaptive function of genetic colour polymorphism in birds. *Biol. Rev. Camb. Philos. Soc.* **79**, 815–848. (doi:10.1017/S1464793104006487)

14. Ducrest A-L, Keller L, Roulin A. 2008 Pleiotropy in the melanocortin system, coloration and behavioural syndromes. *Trends Ecol. Evol.* **23**, 502–510. (doi:10.1016/j.tree.2008.06.001)

15. Galeotti P, Rubolini D, Dunn PO, Fasola M. 2003 Colour polymorphism in birds: causes and functions. *J. Evol. Biol.* **16**, 635–646. (doi:10.1046/j.1420-9101.2003.00569.x)

16. McKinnon JS, Pierotti MER. 2010 Colour polymorphism and correlated characters: genetic mechanisms and evolution. *Mol. Ecol.* **19**, 5101–5125. (doi:10.1111/j.1365-294X.2010.04846.x)

17. Gangoso L, Roulin A, Ducrest AL, Grande JM, Figuerola J. 2015 Morph-specific genetic and environmental variation in innate and acquired immune response in a color polymorphic raptor. *Oecologia* **178**, 1113–1123. (doi:10.1007/s00442- 015-3306-6)

18. Rosenblum EB, Hoekstra HE, Nachman MW. 2004 Adaptive reptile color variation and the evolution of the Mc1r gene. *Evolution* **58**, 1794–1808. (doi:10.1554/03-741)

19. Mullen LM, Hoekstra HE. 2008 Natural selection along an environmental gradient: a classic dine in mouse pigmentation. *Evolution* **62**, 1555–1570. (doi:10.1111/j.1558-5646.2008.00425.x)

20. Rees JL. 2003 Genetics of hair and skin color. *Annu. Rev. Genet.* **37**, 67–90. (doi:10.1146/annurev.genet. 37.110801.143233)

21. Römpler H, Rohland N, Laluzea-Fox C, Willerslev E, Kuznetsova T, Rabeder G, Bertranpetit J, Schöneberg T, Hofreiter M. 2006 Nuclear gene indicates coat-color polymorphism in mammoths. *Science* **313**, 62. (doi:10.1126/science.1128994)

22. Theron E, Hawkins K, Bermingham E, Ricklefs RE, Mundy NI. 2001 The molecular basis of an avian plumage polymorphism in the wild: a melanocortin-1-receptor point mutation is perfectly associated with the melanic plumage morph of the bananaquit, *Coereba flaveola*. *Curr. Biol.* **11**, 550–557. (doi:10.1016/S0960-9822(01)00158-0)

23. Mundy NI, Badcock NS, Hart T, Scribner K, Janssen K, Nadeau NJ. 2004 Conserved genetic basis of a quantitative plumage trait involved in mate choice. *Science* **303**, 1870–1873. (doi:10.1126/science. 1093834)

24. Haas F, Pointer MA, Saino N, Brodin A, Mundy NI, Hansson B. 2009 An analysis of population genetic differentiation and genotype-phenotype association across the hybrid zone of carrion and hooded crows using microsatellites and *MC1R*. *Mol. Ecol.* **18**, 294–305. (doi:10.1111/j.1365-294X. 2008.04017.x)

25. Uy JAC, Moyle RG, Filardi CE, Cheviron ZA. 2009 Difference in plumage color used in species recognition between incipient species is linked to a single amino acid substitution in the melanocortin-1 receptor. *Am. Nat.* **174**, 244–254. (doi:10.1086/600084)

26. Cheviron Z, Hackett SJ, Brumfield RT. 2006 Sequence variation in the coding region of the melanocortin-1 receptor gene (MC1R) is not associated with plumage variation in the blue-crowned manakin (*Lepidoptrix coronata*). *Proc. R. Soc. B* **273**, 1613–1618. (doi:10.1098/rspb.2006.3499)

27. MacDougall-Shackleton EA, Blanchard L, Igdoura SA, Gibbs HL. 2003 Unmelanized plumage patterns in Old World leaf warblers do not correspond to sequence variation at the melanocortin-1 receptor locus (*MC1R*). *Mol. Ecol. Evol.* **20**, 1675–1681. (doi:10.1093/molbev/msg186)

28. Bennett DC, Lamoreux ML. 2003 Pigment gene focus the color loci of mice—a genetic century. *Pigment Cell Res.* **16**, 333–344. (doi:10.1034/j.1600- 0749.2003.00067.x)

29. Kingsley EP, Manceau M, Wiley CD, Hoekstra HE, Borevitz JO. 2009 Melanism in peromyscus is caused by independent mutations in agouti. *PLoS ONE* **4**, e6435. (doi:10.1371/journal.pone.0006435)

30. Poelstra JW, Vijay N, Hoepner MP, Wolf JBW. 2015 Transcriptomics of colour patterning and coloration shifts in crows. *Mol. Ecol.* **24**, 4617–4628. (doi:10.1111/mec.13353)

31. Linnen CR, Kingsley EP, Jensen JD, Hoekstra HE. 2009 On the origin and spread of an adaptive allele in deer mice. *Science* **325**, 1095–1098. (doi:10.1126/ science.1175826)

32. Linnen CR, Poh Y-P, Peterson BK, Barrett RDH, Larson JG, Jensen JD, Hoekstra HE. 2013 Adaptive evolution of multiple traits through multiple mutations at a single gene. *Science* **339**, 1312–1316. (doi:10.1126/science.1233213)

33. Poelstra JW *et al*. 2014 The genomic landscape underlying phenotypic integrity in the face of gene flow in crows. *Science* **344**, 1410–1414. (doi:10.1126/ science.1253226)

34. Ellegren H. *et al*. 2012 The genomic landscape of species divergence in Ficedula flycatchers. *Nature* **491**, 756–760. (doi:10.1038/nature11584)

35. Gill F. 2007 *Ornithology*, 3rd edn. New York, NY: Freeman & Co.

36. Gill FB. 1973 Intra-loral variation in the Mascarene white-eye *Zosterops borbonicus*. *Ornithol. Monogr.* **12**, 1–66. (doi:10.2307/40166689)

37. Bourgeois Y, Bertrand J, Delahaie B, Cornault J, Duval T, Milá B, Thébaud C. 2016 Candidate gene analysis suggests untapped genetic complexity in melanin-based pigmentation in birds. *J. Hered.* **107**, 327–335. (doi:10.1093/jhered/esw017)

38. Bourgeois YXC, Bertrand JAM, Thébaud C, Milá B. 2012 Investigating the role of the melanocortin-1 receptor gene in an extreme case of microgeographical variation in the pattern of melanin-based plumage pigmentation. *PLoS ONE* **7**, e50906. (doi:10.1371/journal.pone.0050906)

39. Bertrand JA, García-Jiménez R, Bourgeois Y, Duval T, Heeb P, Thébaud C, Milá B. 2012 Isolation and characterization of twelve polymorphic microsatellite loci for investigating an extreme case of microgeographical variation in an island bird (*Zosterops borbonicus*). *Conserv. Genet. Resour.* **4**, 323–326. (doi:10.1007/s12686-011-9538-9)

40. Bourgeois YXC *et al*. 2013 Mass production of SNP markers in a nonmodel passerine bird through RAD sequencing and contig mapping to the zebra finch genome. *Mol. Ecol. Resour.* **13**, 899–907. (doi:10.1111/1755-0998.12137)

41. Davey JW, Davey JL, Blaxter ML, Blaxter MW. 2010 RADSeq: next-generation population genetics. *Brief. Funct. Genomics* **9**, 416–423. (doi:10.1093/ bfgp/elq031)
42. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE, Orban L. 2011 A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. *PLoS ONE* **6**, e13979. (doi:10.1371/journal.pone.0019379)

43. Seutin G, White BN, Boag PT. 1991 Preservation of avian blood and tissue samples for DNA analyses. *Can. J. Zool.* **69**, 82–90. (doi:10.1139/z91-013)

44. Griffiths R, Double M, Orr K, Dawson RJG. 1998 A DNA test to sex most birds. *Mol. Ecol.* **7**, 1071–1075. (doi:10.1046/j.1365-294x.1998.00389.x)

45. Cuthill IC, Bennett AFD, Partridge JC, Maier EJ. 1999 Plumage reflectance and the objective assessment of avian sexual dichromatism. *Am. Nat.* **153**, 183–200. (doi:10.1086/303160)

46. Endler JA. 1990 On the measurement and classification of colour in studies of animal colour patterns. *Biol. J. Linn. Soc.* **41**, 315–352. (doi:10.1111/j.1095-8312.1990.tb00839.x)

47. Vorobyev M, Osorio D. 1998 Receptor noise as a determinant of colour thresholds. *Proc. R. Soc. Lond. B* **265**, 351–358. (doi:10.1098/rspb.1998.0302)

48. Cornuault J, Delahaye B, Bertrand JAM, Bourgeois YXC, Milá B, Heeb P, Thébaud C. 2015 Morphological and plumage colour variation in the Réunion grey white-eye (*Aves: Zosterops borbonicus*): assessing the role of selection. *Biol. J. Linn. Soc.* **114**, 459–473. (doi:10.1111/bij.12428)

49. Saino N, Romano M, Rubolini D, Teplitsky C, Ambrosini R, Caprioli M, Canova L, Wakamatsu K. 2013 Sexual dimorphism in melanin pigmentation, feather coloration and its heritability in the barn swallow (*Hirundo rustica*). *PLoS ONE* **8**, e58024. (doi:10.1371/journal.pone.0058024)

50. Ozeki H, Ito S, Wakamatsu K, Thody AJ. 1996 Spectrophotometric characterization of eumelin and phaelelin in hair. *Pigment Cell Res.* **9**, 265–270. (doi:10.1111/j.1600-0749.1996.tb00116.x)

51. Ito S, Nakanishi Y, Valenzuela RK, Brilliant MH, Kolbe L, Wakamatsu K. 2011 Usefulness of alkaline hydrogen peroxide oxidation to analyze eumelin and phaelelin in various tissue samples: application to chemical analysis of human hair melanins. *Pigment Cell Melanoma Res.* **24**, 605–613. (doi:10.1111/j.1755-148X.2011.00864.x)

52. Wakamatsu K, Ito S, Rees JL. 2002 The usefulness of 4-amino-3-hydroxyphenylalanine as a specific marker of phaelelin. *Pigment Cell Res.* **15**, 225–232. (doi:10.1034/j.1600-0749.2002.02009.x)

53. Bertrand JAM *et al.* 2014 Extremely reduced dispersal and gene flow in an island bird. *Heredity* **112**, 190–196. (doi:10.1038/hdy.2013.91)

54. Bertrand JAM, Delahaye B, Bourgeois YXC, Duval T, García-Jiménez R, Cornuault J, Pujol B. 2016 The role of selection and historical factors in driving population differentiation along an elevational gradient in an island bird. *J. Evol. Biol.* **29**, 824–836. (doi:10.1111/jeb.12829)

55. Hadfield JD, Richardson DS, Burke T. 2006 Towards unbiased parentage assignment: combining genetic, behavioural and spatial data in a Bayesian framework. *Mol. Ecol.* **15**, 3715–3730. (doi:10.1111/j.1365-294X.2006.03050.x)

56. Cornetti L, Valente LM, Dunning LT, Quan X, Black RA, Hébert OH, Savolainen V. 2015 The genome of the 'great speciator' provides insights into bird diversification. *Genome Biol. Evol.* **7**, 2680–2691. (doi:10.1093/gbe/evv168)

57. Li H, Durbin R. 2009 Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 1754–1760. (doi:10.1093/bioinformatics/btp324)

58. Schwartz S, Kent W, Smit A. 2003 Human–mouse alignments with BLASTZ. *Genome Res.* **13**, 103–107. (doi:10.1101/gr.809403)

59. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer M, Marth G, Abecasis G, Durbin R. 2009 The sequence alignment/map format and SAMtools. *Bioinformatics* **25**, 2078–2079. (doi:10.1093/bioinformatics/btp352)

60. Kofler R, Pandey RV, Schlötterer C. 2011 PoPoolation2: Identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq). *Bioinformatics* **27**, 3435–3436. (doi:10.1093/bioinformatics/btr589)

61. Patterson N, Price AL, Reich D. 2006 Population structure and eigenanalysis. *PLoS Genet.* **2**, e190. (doi:10.1371/journal.pgen.0020190)

62. Huber W *et al.* 2015 Orchestrating high-throughput genomic analysis with Conductor. *Nat. Methods* **12**, 115–121. (doi:10.1038/nmeth.3252)

63. Excoffier L, Lischer HEL. 2010 Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. *Mol. Ecol. Resour.* **10**, 564–567. (doi:10.1111/j.1755-0998.2010.02847.x)

64. Gautier M. 2015 Genome-wide scan for adaptive divergence and association with population-specific covariates. *Genetics* **201**, 1555–1579. (doi:10.1534/genetics.115.181453)

65. Purcell S *et al.* 2007 PLINK: a tool set for whole-genome association and population-based linkage analyses. *Am. J. Hum. Genet.* **81**, 559–575. (doi:10.1086/519759)

66. Manichaikul A, Mychaleckyj JC, Rich SS, Daly K, Sale M, Chen W-M. 2010 Robust relationship inference in genome-wide association studies. *Bioinformatics* **26**, 2867–2873. (doi:10.1093/bioinformatics/btq559)

67. Li M, Boehnke M, Abecasis GR. 2006 Efficient study designs for test of genetic association using sibship data and unrelated cases and controls. *Am. J. Hum. Genet.* **78**, 778–792. (doi:10.1086/503711)

68. Doran AG, Creevey CJ. 2013 Snpatd: easy and rapid annotation of results from de novo snp discovery projects for model and non-model organisms. *BMC Bioinform.* **14**, 45. (doi:10.1186/1471-2105-14-45)

69. Rasmussen MD, Hubisz MJ, Gronau I, Siepel A, Coop G. 2014 Genome-wide inference of ancestral recombination graphs. *PLoS Genet.* **10**, e1004342. (doi:10.1371/journal.pgen.1004342)

70. Prum RO, Dyck J. 2003 A hierarchical model of plumage: morphology, development, and evolution. *J. Exp. Zool. B. Mol. Dev. Evol.* **298B**, 73–90. (doi:10.1002/jez.b.27)

71. Weber BHF *et al.* 2002 Inactivation of the murine X-linked juvenile retinoschisis gene, Rs1 h, suggests a role of retinoschisis in retinal cell layer organization and synaptic structure. *Proc. Natl Acad. Sci. USA* **99**, 6222–6227. (doi:10.1073/pnas.092528599)

72. Delevoye C *et al.* 2009 AP-1 and KIF13A coordinate endosomal sorting and positioning during melanosome biogenesis. *J. Cell Biol.* **187**, 247–264. (doi:10.1083/jcbc.200907122)

73. Theos A, Tenza D, Martina J. 2005 Functions of adaptor protein (AP)-3 and AP-1 in tyrosinase sorting from endosomes to melanosomes. *Mol. Biol. Cell* **16**, 5356–5372. (doi:10.1091/mbc.E05)

74. Feng L *et al.* 1999 The beta3a subunit gene (Ap3b1) of the AP-3 adaptor complex is altered in the mouse hypopigmentation mutant pearl, a model for Hermansky-Pudlak syndrome and night blindness. *Hum. Mol. Genet.* **8**, 323–330. (doi:10.1093/hmg/8.2.323)

75. Hoek KS *et al.* 2008 Novel MITF targets identified using a two-step DNA microarray strategy. *Pigment Cell Melanoma Res.* **21**, 665–676. (doi:10.1111/j.1755-148X.2008.00505.x)

76. Gwynn B, Smith RS, Rowe LB, Taylor BA, Peters LL. 2006 A mouse TRAPP-related protein is involved in pigmentation. *Genomics* **88**, 196–203. (doi:10.1016/j.ygeno.2006.04.002)

77. Ohbayashi N, Fukuda M. 2012 Role of Rab family GTPases and their effectors in melanosomal logistics. *J. Biochem.* **151**, 343–351. (doi:10.1093/jmb/mys009)

78. Kloer DP *et al.* 2010 Assembly of the biogenesis of lysosome-related organelles complex-3 (BLOC-3) and its interaction with Rab9. *J. Biol. Chem.* **285**, 7794–7804. (doi:10.1074/jbc.M109.069088)

79. Fairbank PD, Lee C, Ellis A, Hildebrand JD, Gross JM, Wallingford JB. 2006 Shroom2 (APXL) regulates melanosome biogenesis and localization in the retinal pigment epithelium. *Development* **133**, 4109–4118. (doi:10.1242/dev.02563)

80. Palmisano I *et al.* 2008 The ocular albinism type 1 protein, an intracellular G protein-coupled receptor, regulates melanosome transport in pigment cells. *Hum. Mol. Genet.* **17**, 3487–3501. (doi:10.1093/hmg/ddn242)

81. Incerti B *et al.* 2000 Oa1 knock-out: new insights on the pathogenesis of ocular albinism type 1. *Hum. Mol. Genet.* **9**, 2781–2788. (doi:10.1093/hmg/9.19.2781)

82. Uy JAC *et al.* 2016 Mutations in different pigmentation genes are associated with parallel melanism in island flycatchers. *Proc. R. Soc. B* **283**, 2115–2118. (doi:10.1098/rspb.2016.0731)

83. Roulin A, Eremasi G, Bize P, Gasparini J, Piault R, Ducreux A-L. 2011 Pale and dark reddish melanic tawny owls differentially regulate the level of blood circulating POMC prohormone in relation to environmental conditions. *Oecologia* **166**, 913–921. (doi:10.1007/s00442-011-1955-7)

84. Flint J, Mackay TFC. 2009 Genetic architecture of quantitative traits in mice, flies, and humans. *Genome Res.* **19**, 723–733. (doi:10.1101/gr.086660.108.19)

85. Joron M *et al.* 2011 Chromosomal rearrangements maintain a polymorphic supergene controlling butterfly mimicry. *Nature* **477**, 203–206. (doi:10.1038/nature10341)

86. Mallet J, Barton N. 1989 Strong natural selection in a warning color hybrid zone. *Evolution* **43**, 421–431. (doi:10.2307/2409217)

87. Turner JRG. 1977 Butterfly mimicry: the genetic evolution of an adaptation. *Evol. Biol.* **11**, 163–206. (doi:10.1007/978-1-4615-6953-4_3)

88. Bai J, Xie X, Lei Y, An G, He L, Lv X. 2014 Ocular albinism type 1 induced melanoma cell migration is mediated through the RAS/RAF/MEK/ERK signaling pathway.
pathway. Mol. Med. Rep. 10, 491-495. (doi:10.3892/mmrr.2014.2154)

89. Drabek K, van de Peppel J, Eijken M, van Leeuwen JPTM. 2011 GPM6B regulates osteoblast function and induction of mineralization by controlling cytoskeleton and matrix vesicle release. *J. Bone Miner. Res.* **26**, 2045–2051. (doi:10.1002/jbmr.435)

90. Fjorback AW, Müller HK, Wiborg O. 2009 Membrane glycoprotein M6B interacts with the human serotonin transporter. *J. Mol. Neurosci.* **37**, 191–200. (doi:10.1007/s12031-008-9092-4)

91. Sinervo B, Svensson E. 2002 Correlational selection and the evolution of genomic architecture. *Heredity* **89**, 329–338. (doi:10.1038/sj.hdy.6800148)

92. Abzhanov A, Extavour CG, Groover A, Hodges SA, Hoekstra HE, Kramer EM, Monteiro A. 2008 Are we there yet? Tracking the development of new model systems. *Trends Genet.* **24**, 353–360. (doi:10.1016/j.tig.2008.04.002)

93. Teshima KM, Przeworski M. 2006 Directional positive selection on an allele

of arbitrary dominance. *Genetics* **172**,
713–718. (doi:10.1534/genetics.105.
044065)

94. Charlesworth B, Charlesworth D. 2010 Elements of evolutionary genetics. *Bioscience* **61**, 409–411. (doi:10.1525/bio.2011.61.5.12)

95. Bourgeois YXC et al. 2017 Data from: A novel locus on chromosome 1 underlies the evolution of a melanic plumage polymorphism in a wild songbird. Dryad Digital Repository. (http://dx.doi.org/10.5061/dryad.1q9g7)