Genomic and plumage variation in Vermivora hybrids

Abstract

Hybrids with different combinations of traits can be used to identify genomic regions that underlie phenotypic characters important to species identity and recognition. Here, we explore links between genomic and plumage variation in Blue-winged Warbler x Golden-winged Warbler (Vermivora cyanoptera x V. chrysoptera) hybrids, which have traditionally been categorized into 2 discrete types. “Lawrence’s” hybrids are yellow overall, similar to Blue-winged Warblers, but exhibit the black throat patch and face mask of Golden-winged Warblers. “Brewster’s” hybrids are similar to Golden-winged Warblers, but lack the black throat patch and face mask, and sometimes have yellow on their underparts. Previous studies hypothesized that (1) first generation hybrids are of the Brewster’s type and can be distinguished by the amount of yellow on their underparts, and that (2) the throat patch/mask phenotype is consistent with Mendelian inheritance and controlled by variation in a locus near the Agouti-signaling protein (ASIP) gene. We addressed these hypotheses using whole genome re-sequencing of parental and hybrid individuals. We found that Brewster’s hybrids had genomic hybrid index scores indicating this phenotype can arise by majority ancestry from either parental species, that their plumage varied in levels of carotenoid pigmentation, and individuals captured in multiple years grew consistently less yellow over time. Variation in carotenoid pigmentation showed little relationship with genomic hybrid index score and is thus inconsistent with previous hypotheses that first generation hybrids can be distinguished by the amount of yellow in their plumage. Our results also confirm that variation near ASIP underlies the throat patch phenotype, which we refined to an ~10–15 Kb region upstream of the coding sequence. Overall, our results support the notion that traditional categorization of hybrids as either Lawrence’s or Brewster’s oversimplifies continuous variation in carotenoid pigmentation, and its inferred underlying genetic basis, and is based primarily on one discrete trait, which is the throat patch/mask phenotype.

Resumen

Los híbridos con diferentes combinaciones de rasgos pueden ser usados para identificar regiones genómicas que subyacen a los caracteres fenotípicos importantes para la identidad y el reconocimiento de las especies. Aquí exploramos los vínculos entre la variación genómica y el plumaje en híbridos de Vermivora cyanoptera x V. chrysoptera, que han sido tradicionalmente categorizados en dos tipos discretos. Los híbridos “Lawrence” son muy amarillos en general, similares a V. cyanoptera, pero exhiben el parche gular y la mascarilla negros de V. chrysoptera. Los híbridos de “Brewster” son similares a V. chrysoptera, pero no tienen el parche gular y la mascarilla negros y a veces presentan amarillo en las partes inferiores. Estudios previos hipotetizaron que (1) los híbridos de la primera generación son del tipo Brewster y pueden ser distinguidos por la cantidad de amarillo en sus partes inferiores, y que (2) el fenotipo del parche gular y la mascarilla es consistente con la herencia Mendeliana y es controlado por la variación en un locus cerca del gen de la proteína de señalización de Agutí (ASIP, por sus siglas en inglés). Abordamos estas hipótesis usando re-secuenciación genómica completa de individuos parentales e híbridos. Encontramos que los híbridos de Brewster tuvieron puntajes del índice genómico híbrido que indicaron que este fenotipo puede surgir por una mayoría ancestral a partir de cualquiera de las especies parentales, que sus plumajes variaron en los niveles de pigmentación de carotenoides, y que los individuos capturados en múltiples años desarrollaron consistentemente menos amarillo a lo largo del tiempo. La variación en la pigmentación de carotenoides mostró poca relación con el índice genómico híbrido y es por ende consistente con hipótesis previas que establecen que los híbridos de la primera generación pueden ser distinguidos por la cantidad de amarillo en sus plumajes. Nuestros resultados también confirman que la variación cerca de la ASIP subyace el fenotipo de parche gular, la cual nosotros refinamos a una región aproximada de 10–15 Kb hacia arriba de la secuencia de codificación. En general, nuestros resultados apoyan la noción de que la categorización tradicional de los híbridos como Lawrence o Brewster es una simplificación exagerada de la variación continua en la pigmentación de carotenoides y de su base genética subyacente inferida, y está basada principalmente en un solo rasgo discreto, que es el fenotipo de parche gular/mascarilla.

Among birds, wood warblers are among the most variable in their vibrant plumage pigmentation. The evolutionary and genetic mechanisms that underlie plumage differences are not well known. We address this here by sequencing natural hybrids between Blue-winged and Golden-winged warblers.

INTRODUCTION

Naturalists have long been fascinated by the incredible diversity in avian plumage coloration (Cuthill et al. 2017, Funk and Taylor 2019), which can vary dramatically even between closely related species (e.g., Poelstra et al. 2014). Plumage coloration is hypothesized to play a strong role in many facets of avian evolution, including sexual selection, species recognition, feather durability, thermoregulation, and crypsis (Burtt 1986, Savalli 1995). Variation in feather coloration can result from the physical structure of the feather (e.g., McCoy et al. 2018) and/or the pigment it contains (e.g., Poelstra et al. 2014), which may be influenced by an individual’s diet (e.g., Hill 1992). Thus, understanding the evolutionary, genetic, and developmental mechanisms that produce and maintain species differences in plumage coloration is of broad interest.

In closely related avian groups that differ strongly in coloration, it is often difficult to identify functional genetic variation underlying divergent plumages, because the independent accumulation of many mutations across the genome during species divergence masks the few loci that are responsible for these phenotypic differences. Systems with substantial natural hybridization, where distinct lineages interbreed and produce admixed offspring (hybrids), can be used to address this problem. Given multiple generations of admixture and backcrossing, hybrid chromosomes are recombinant tracts of ancestry from the parental species. Thus, the power to detect associations between discrete phenotypes and their underlying genetic variants is increased in studies of hybrid swarms and other taxa with low differentiation and/or substantial ongoing admixture (Harrison 1990, Rieseberg and Buerkle 2002, Brelsford et al. 2017, Campagna et al. 2017).

The Blue-winged Warbler x Golden-winged Warbler (Vermivora cyanoptera x V. chrysoptera) mosaic hybrid zone system is a highly tractable model for genetic studies of plumage traits. The parental species are sister taxa (Lovette et al. 2010) with striking differences in feather pigmentation patterns (Figure 1B). However, despite these clear phenotypic differences, genome-wide differentiation between Blue-winged and Golden-winged Warblers is low (Toews et al. 2016), likely because they have long hybridized in regions where their breeding ranges overlap in eastern North America.

FIGURE 1.

(A) Genomic admixture and plumage yellow scores for hybrids sequenced in this study. Hybrid index is an estimate of the proportion of the genome inherited from the Golden-winged Warbler and interspecific heterozygosity is the proportion of the genome heterozygous with an allele inherited from both parental species. Golden-winged Warblers have a hybrid index score of 1 and Blue-winged Warblers have a hybrid index of 0. Double circles denote individuals that were captured and scored for plumage in multiple years (see Results). For these individuals, we plotted the average plumage score. (B) Plumage types for parental species and hybrids. Symbols for hybrids and hybrid index scores for parental species are consistent with their representation in panel (A). Warbler illustrations by Liz Clayton Fuller.

The parental species can be readily distinguished by 3 main plumage characters (Parkes 1951): (1) presence/absence of the black throat patch and face mask, which are correlated (i.e. individuals with a black throat also have a black mask, and individuals with a plain throat lack the black mask); (2) body color, especially the amount of yellow (carotenoid-based pigments) on the underparts; and (3) the width, separation, and color of the wingbars. Although Vermivora hybrids display a range of intermediate phenotypes, 2 forms are common and discrete enough from the parental types to have been named (Figure 1B). “Lawrence’s” hybrids are similar to Blue-winged Warblers (i.e. yellow overall, with 2 narrow white wing bars) but have the black throat patch and face mask, similar to Golden-winged Warblers. “Brewster’s” hybrids, by contrast, lack a black throat patch, have little to no yellow on the underparts, and commonly have partially separated yellowish wing bars. Variation within each type has been observed, particularly in the amount of yellow on the underparts among Brewster’s hybrids (Parkes 1951, Short 1963). Other hybrid configurations resembling Golden-winged Warblers with varying amounts of yellow on the breast (typically referred to as “backcrossed hybrids” or “intergrades”) also occur.

Nichols (1908) speculated about the genetics of these plumage characters in Vermivora warblers and hypothesized that they are inherited in a Mendelian fashion, with plain throat (pp) and absence of white (i.e. yellow) underparts (ww) being homozygous recessive. First generation (F1) hybrids would then be of the Brewster’s type, with a plain throat and white underparts (WwPp). Taking into account variation in the amount of yellow on the underparts of Brewster’s individuals observed in museum collections, Parkes (1951) instead speculated that the allele for white underparts (W) is actually incompletely dominant, and heterozygous F1 Brewster’s warblers would have partially yellow underparts. If the white underparts allele is incompletely dominant, then Brewster’s individuals with all-white underparts would have to be homozygous for the white allele (WW), likely a result of backcrossing with Golden-winged Warblers. Lawrence’s warblers (with yellow underparts and black throat) would then have to be homozygous for both alleles (wwpp), which could arise in F2 and later generations. Parkes (1951) also noted that abundances of the Brewster’s and Lawrence’s types in zones of contact seem to be consistent with the expected frequencies of hybrid types upon backcrossing of F1 Brewster’s individuals.

Toews et al. (2016) found that the genomes of the parental species are nearly homogenous except for 6 highly differentiated regions, 4 of which are near genes previously linked to plumage traits. In particular, genotyping of a variant near the Agouti-signaling protein (ASIP) gene revealed a perfect correlation with presence/absence of the black throat patch in a large sample of individuals including hybrids. Consistent with the predictions of Nichols (1908), ASIP appeared to be inherited and expressed as a Mendelian trait, with the black throat recessive to plain throat color in heterozygotes. This variant is located in a moderately large (>50 Kb) divergent region of the genome upstream of the ASIP gene and its homologous chicken promoter regions, which presumably contains the causal promoter in Vermivora warblers. Variants near this gene have been linked to feather pigmentation in other avian systems (capuchino seedeaters: Campagna et al. 2017, munias: Stryjewski and Sorenson 2017), where divergence peaks were similarly large (~90–100 Kb), indicating that ASIP likely plays a role in plumage divergence, but the causal mutations within divergent regions that underlie species differences are difficult to identify.

In this study, we employed whole genome sequencing to investigate links between phenotypic and genomic variation in hybrid Vermivora individuals. We first directly tested whether Brewster’s hybrid genomes are consistent with Parkes’ (1951) plumage hypothesis: Do early generation Vermivora hybrids express more carotenoid-based pigments, and white-breasted hybrids result via backcrossing with Golden-winged Warblers? Or, do these early generation hybrids instead express few carotenoids, as predicted by Nichols (1908)? Next, we took advantage of variation and recombination in the genomes of phenotypic hybrids to refine the genomic region containing the ASIP promoter controlling the throat patch phenotype. Previous studies using whole genome sequencing only focused on parental individuals, which did not allow for this type of refined analysis.

METHODS

Sampling and Sequencing

We obtained blood samples, taken from the brachial vein, from breeding, territorial male warblers captured from 2015 to 2017 at 14 localities in 4 states spanning the hybrid zone between Vermivora warblers (Appendix Table 2). Blood samples were stored in a 0.5% sodium dodecyl sulfate (SDS) buffer in the field and stored at room temperature prior to DNA extraction. Samples spanned a range of different parental and hybrid phenotypes. We sampled 20 individuals at the ends of the plumage spectrum (i.e. parental individuals; n = 10 Blue-winged Warblers, n = 10 Golden-winged Warblers). We also sampled 9 phenotypic hybrids: 6 classified as the Brewster’s warbler hybrid phenotype, 2 classified as “intergrades”, and 1 Lawrence’s warbler hybrid phenotype (Figure 1). Samples from the different phenotypic classes were expressly chosen to represent individuals in that given class, except for Brewster’s hybrids, which were chosen arbitrarily from the available samples. Data for these hybrid samples have been published along with 2 previous studies (Toews et al. 2016, 2019), although the hybrid genomes have not been explicitly analyzed. Sequencing data for the one Lawrence’s hybrid individual (Table 1) has not been previously published.

TABLE 1.

Sample information for individuals sequenced in this study. Individual sampling localities and plumage scores can be found in Appendix Tables 2 and 3.

Plumage phenotype	Record number	Plumage yellow score (out of 31)	Accession numbers
Golden-winged Warbler a	Multiple (n = 10)	0–2	SAMN05223-489, -493, 497, -499, -507, -608, -512, -514, -518, -526
Blue-winged Warbler a	Multiple (n = 10)	29–31	SAMN05223-487, -492, -495, -501, -503, -506, -509, -519, -523, -528
Brewster’s b	MIH02	9, 18	SAMN12573390
Brewster’s b	PAH01	4, 12	SAMN12573391
Brewster’s	PAH06	10	SAMN12573392
Brewster’s b	WIH03	4, 5, 6	SAMN12573393
Brewster’s a	215081706	10	SAMN05223527
Brewster’s a	215081713	10	SAMN05223522
Intergrade a	215081718	14	SAMN05223502
Intergrade a	215081729	7	SAMN05223533
Lawrence’s	PE25T03	29	SAMN14431356

^a Library preparation using PCRFree kit.

^b Individuals captured in multiple years; plumage yellow scores are reported for each year (MIH02 and PAH01 were captured in 2015 and 2016, WIH03 was captured in 2015, 2016, and 2017).

For the plumage scores, we followed the scoring criteria developed by (Gill 1980), with the addition of 3 characters: mask, eye-line, and mustache color (Toews et al. 2017). Birds were scored for the 11 characters in the hand (the parental taxa) as well as from digital photos (the hybrids), and all scores were quantified by a single researcher (D.P.L.T.) Photographs are available for all samples in the Dryad dataset at https://doi.org/10.5061/dryad.tx95x69tk. Because we used the plumage scores to analyze yellow coloration, we report the modified “sum yellow” plumage score, implemented by Toews et al. (2017), which only includes carotenoid-based plumage characters, and inverted it by subtracting from the maximum value (31) such that higher scores represent more yellow plumage (Appendix Table 3).

We extracted DNA from all samples using UPrep spin columns (Genesee Scientific, Rochester, New York, USA), following the Qiagen (Germantown, Maryland, USA) buffer blood extraction protocols. We combined data from different sequencing projects that utilized 2 library preparation methods, which differ slightly. For 24 samples, we used the Illumina (San Diego, California, USA) TruSeq PCRFree kit, following the protocols to generate libraries with a 350 base pair (bp) insert. For the remaining 9 individuals, DNA concentrations were lower and, therefore, we used the Illumina TruSeq Nano kit—which includes an 8-cycle polymerase chain reaction (PCR) enrichment—also selecting 350 bp insert sizes. We sequenced these libraries across separate lanes of an Illumina NextSeq using the paired-end 150 bp sequencing chemistry. In some cases, we combined Vermivora samples with other wood warblers from other projects, but consistently included 24 individuals per sequencing lane in order to generate comparable coverage across samples.

Sequence Processing and Alignments

We used AdapterRemoval 2.1.7 to trim Ns from 5’ and 3’ ends and to collapse overlapping read pairs. Reads shorter than 20 bp after trimming were discarded. We aligned both un-collapsed and collapsed reads to the Yellow-rumped Warbler (Setophaga coronata) genome assembly (Toews et al. 2016) using the “very sensitive local” setting in Bowtie 2 2.3.5.1 (Langmead and Salzberg 2012). We marked PCR duplicates using default settings with Picard MarkDuplicates 2.20.4 (http://broadinstitute.github.io/picard/). Because we were interested in genomic variation underlying plumage differences between the parental species, we restricted all analyses to the 6 previously described genomic regions of high differentiation (median size ~30 Kb) described in Toews et al. (2016) (see Appendix Table 4).

Genotype Likelihoods and Genotyping

Due to low coverage in our resequencing data, we employed a method to estimate genotype likelihoods from our alignments, rather than calling genotypes directly, which can inflate uncertainty by calling genotypes from sequencing errors. To do this, we used the Genome Analysis Toolkit (GATK) model in ANGSD 0.929 (Korneliussen et al. 2014) to infer major and minor allele frequencies and removed sites with a P-value > 1e–6.

Some of the tools we used in our analyses, however, required called genotypes as input. Thus, we also called genotypes from our alignments using GATK’s HaplotypeCaller 3.8 (DePristo et al. 2011) with default settings and retained sites with a minimum depth of 4 using vcftools 0.1.17 (Danecek et al. 2011).

Genomic Variation in Highly Divergent Regions

We used principal component analysis to visualize genomic variation in the 6 highly divergent regions. To do this, we used the genotype likelihoods to estimate covariance matrices with PCAngsd 0.981 (Meisner and Albrechtsen 2018) using default settings. We then used the eigen function in base R 3.5.1 to compute and plot the first 2 eigenvectors.

Hybrid Admixture Estimation

We used the R package introgress (Gompert and Buerkle 2010) on the genotype data to calculate 2 estimates of genomic admixture for hybrids: hybrid index and interspecific heterozygosity. Hybrid index is the proportion of the genome (in this case a subset; see below) inherited from one parental species, in this case phenotypic Golden-winged Warblers. Thus, F1 hybrids have a hybrid index score of 0.5, hybrids backcrossed into Golden-winged Warblers will have a hybrid index >0.5, and hybrids backcrossed into Blue-winged Warblers will have a hybrid index <0.5. Interspecific heterozygosity is the proportion of the genetic markers used in the hybrid index score that are heterozygous, with an allele from each parental species (i.e. Blue-winged Warbler/Golden-winged Warbler), where for F1 hybrids interspecific heterozygosity is equal to one.

It is important to note that since we necessarily restricted our genotyping to the 6 small divergent regions (introgress requires nearly fixed markers), these “hybrid index” scores are not representative of admixture across the whole genome. This contrasts with the convention in other hybrid zone studies that exhibit more genome-wide differentiation (e.g., Trier et al. 2014, Toews et al. 2018), where the hybrid index is a better summary of genome-wide admixture patterns. Instead, in Vermivora, the scores are meant to represent the expected genomic variation underlying species differences in known plumage traits, and are used to roughly distinguish the ancestry of the different hybrid classes (Toews et al. 2016). For this analysis, we used vcftools to calculate site-specific fixation indices (F_ST) and retained only sites in divergent regions that were highly differentiated (F_ST > 0.7) between phenotypically Blue-winged Warbler and Golden-winged Warbler individuals.

Refining the ASIP Causal Region

We used 2 methods to further resolve the genomic region underlying the throat patch phenotype. Our first method assayed variation in nucleotide diversity across the warbler scaffold containing ASIP (warbler scaffold 299) and relied on the following assumptions: (1) the presence of the black throat patch is homozygous recessive, (2) the recessive allele is fixed in Golden-winged Warblers and the dominant allele is fixed or at high frequency in Blue-winged Warblers, and (3) Brewster’s hybrids are heterozygous at the causal site. Our first assumption implies that Golden-winged Warblers, intergrades, and Lawrence’s hybrids are homozygous recessive for the causal allele since these individuals display the black throat patch, which is consistent with previous genomic analysis (Toews et al. 2016) and studies of hybrid plumage (Parkes 1951). Our second assumption implies that some Blue-winged Warblers may be heterozygous for the black throat allele although they exhibit a plain throat. Assuming the parental species exhibit fixed differences in the throat patch allele is consistent with range-wide sampling and genotyping of a marker in the candidate region (Toews et al. 2016). Our last assumption is also consistent with Toews et al. (2016) and is conservative, since backcrossed hybrids with Brewster’s plumage could be homozygous for the dominant plain throat allele.

Following these assumptions, we used vcftools to calculate individual nucleotide diversity (π) for sites along the scaffold containing ASIP (warbler scaffold 299). Since we employed this as a within-individual measure, heterozygous sites will have π = 1 and homozygous sites will have π = 0. We then retained sites where all genotyped Brewster’s individuals were heterozygous and the remaining individuals (Blue-winged Warblers, Golden-winged Warblers, Lawrence’s, intergrades) were homozygous, with the expectation that these sites will cluster around the variant underlying the throat patch phenotype (i.e. there will be increased diversity in Brewster’s individuals). We computed and plotted kernel density estimates for these sites using the density function in R.

Our second method included finding possible recombination breakpoints in hybrids using multivariate analysis. In this case, the first eigenvector in the principal component analysis of the ASIP-containing scaffold strongly separates Blue-winged Warblers from Golden-winged Warblers. We therefore conducted principal component analyses iteratively along this scaffold to refine the causal region. In particular, we focused on the first eigenvector score of the Lawrence’s warbler (the hybrid exhibiting mostly Blue-winged Warbler-like plumage traits that also has the black throat patch of Golden-winged Warblers). If the throat patch allele is located within this scaffold, we expect the Lawrence’s warbler to cluster with Blue-winged Warblers (plain throat) outside the causal region, but to only cluster with Golden-winged Warblers (black throat) within the causal region. For this analysis, we computed covariance matrices in 5 Kb windows within the highly divergent region in the scaffold containing ASIP (warbler scaffold 299) using the methods described above.

RESULTS

Sequencing and Genotyping Rates

After trimming, we retained an average of ~32.5 million reads per individual (range: ~20–47 million reads). Mean coverage per individual ranged from ~4 to 11X across the scaffolds we used in our analysis. Using our genotype calling approach, we obtained 1,434,494 single nucleotide polymorphisms (SNPs) within the 6 divergent genomic regions. After filtering, we used 960 SNPs in our hybrid admixture analyses (only SNPs in divergent windows with depth > 4 and F_ST > 0.7). The genotyping rate for these 960 SNPs was high across individuals (mean = 93%). We used 30,747 SNPs in our nucleotide diversity analyses of warbler scaffold 299, which contains ASIP (mean genotyping rate = 90%). Using our genotype likelihood-based approach, we retained between 586 and 23,077 sites per scaffold to estimate covariance matrices within divergent regions, and between 394 and 627 sites within the 5 Kb windows along scaffold 299 (Appendix Table 4).

Vermivora Hybrids are Highly Variable

We observed considerable variation within divergent genomic regions and in phenotypes of hybrid individuals. Our estimates of hybrid index ranged from 0.11 (Lawrence’s) to 0.95 (Brewster’s), and interspecific heterozygosity ranged from 0.11 (Brewster’s) to 0.58 (Brewster’s) (Figure 1A). Although some hybrids had intermediate hybrid index scores, it is likely that none are F1 hybrids since heterozygosity was always lower than 0.6, even among the Brewster’s hybrids.

We found admixture scores and interspecific heterozygosity for hybrid individuals in the present study consistent with those same individuals genotyped previously by Toews et al. (2016) in a low-resolution restriction fragment length polymorphism (RFLP) assay (overlap of n = 4 individuals between studies; Appendix Figure 5). Thus, despite the low coverage of our sequencing data, our ability to detect F1s is not likely hindered by an inability to call heterozygous genotypes, but rather the fact that there were few F1 hybrids in our sample. Toews et al. (2016) identified one individual that was heterozygous across all 6 RFLP sites. This individual had a Brewster’s phenotype and may represent an F1 hybrid (see discussion below), but was not included in the present study as we did not want to bias our sampling of Brewster’s warblers based on previous results.

Consistent with Lawrence’s hybrids resulting from backcrossing into Blue-winged Warblers, the Lawrence’s individual had a low hybrid index. Similarly, the intergrades had high hybrid indices, consistent with backcrossing into Golden-winged Warblers. Hybrid index scores for Brewster’s hybrids were quite variable (range: 0.34–0.95), suggesting the Brewster’s phenotype can result as a consequence of backcrossing into either parental species.

Hybrid plumage scores ranged from 26 (Lawrence’s) to 1 (intergrade). Based on visual assessment of the data, plumage yellow score does not appear to be related to genomic hybrid index score (Table 1, Figure 1A). We also noted an unpredicted pattern in 2 hybrids that were captured and photographed in 2 consecutive years, and in one that was captured and photographed in 3 consecutive years (all Brewster’s hybrids). Surprisingly, for each of these hybrids, plumage scores changed each year—in some cases dramatically—as a result of the loss of yellow coloration (Figure 2). We are unaware of other documented instances of such significant inter-annual changes in adult breeding plumage in these or other wood warblers.

FIGURE 2.

Unexpected variation in carotenoid-based plumage for hybrids captured in multiple years. (A) Plumage yellow scores, (B) a Brewster’s hybrid (MIH02, hybrid index = 0.5, interspecific heterozygosity = 0.29) photographed in 2015, and (C) in 2016. Individual IDs reflect sampling locality: MI = Michigan, PA = Pennsylvania, WI = Wisconsin.

Principle component analyses of divergent regions were also consistent with variability among hybrids (Figure 3). For example, in some regions, Brewster’s hybrids clustered intermediately between parental types (e.g., warbler scaffold 120), while in other regions, some Brewster’s hybrids clustered intermediately and some clustered with either parental type (e.g., warbler scaffold 24). The intergrades usually clustered with Golden-winged Warblers, although there are some exceptions (i.e. warbler scaffolds 120, 24, 653). The Lawrence’s hybrid clustered consistently with Blue-winged Warblers, except for along the region of the scaffold containing ASIP (warbler scaffold 299). Here, its eigenvector score was closer to Golden-winged Warblers, consistent with this region underlying the throat patch phenotype.

FIGURE 3.

Results of principle components analysis of genomic variation within the 6 divergent regions for parental types (BWWA = Blue-winged Warblers, GWWA = Golden-winged Warblers) and hybrids (see Appendix Table 4 for inclusive regions). Warbler scaffold 299 contains the ASIP gene.

ASIP Promotor Region Narrowed to ~10–15 Kb

Our nucleotide diversity analysis of warbler scaffold 299 revealed that sites where Brewster’s individuals were heterozygous, and the other plumage types (Lawrence’s, intergrades, Blue-winged Warbler, Golden-winged Warbler) were homozygous, clustered at the far end of warbler scaffold 299, near position 390,000 (Figure 4B). Our iterative principle component analysis was also consistent with this finding. In 5 Kb windows before position 385,000 and after position 400,000, the Lawrence’s hybrid clusters intermediate to parental types. However, within the windows between 385,000 and 400,000, the Lawrence’s hybrid clusters with Golden-winged Warblers and the intergrades (Figure 4A), suggesting this ~10–15 Kb region contains the sequence that underlies the throat patch phenotype in Vermivora warblers.

FIGURE 4.

Genomic variation on the warbler scaffold containing the ASIP gene. (A) Results of our iterative principle components analysis of 5 Kb windows within the presumed causal throat patch region, (B) density plot of sites along this scaffold where Brewster’s hybrids were heterozygous (π = 1) and other plumage types were homozygous (π = 0). This warbler scaffold aligns to Zebra Finch chromosome 20:1,483,796–1,967,099 (Toews et al. 2016).

DISCUSSION

Here, we used whole genome re-sequencing to explore links between genomic and phenotypic variation in Vermivora hybrids. The variation we observed suggests that the genomic underpinnings of the hybrid plumage types in this system are not straightforward, despite divergence being restricted to a few genomic regions in the parental species. One exception is for the region controlling the throat patch phenotype, which we have narrowed to an ~10–15 Kb window upstream of the ASIP gene.

Genomic Variation and Discrete Hybrid Types

Carotenoid-based plumage scores were variable among hybrid types (Brewster’s, Lawrence’s, intergrades) and did not have a strong relationship with the degree of genetic admixture (Figure 1A). This is in agreement with the continuous variation in hybrid color characters observed by Short (1963, 1969), who argued that it is not appropriate to categorize Vermivora hybrids into discrete classes (i.e. Lawrence’s or Brewster’s). Our results add additional nuance to this argument by demonstrating that genomic variation in hybrids of the Brewster’s type is continuous and that phenotypic and genomic data support only one discrete trait for hybrid types in this system, which is the throat patch phenotype. Our results suggest that (1) the Lawrence’s individual we sequenced is a multigenerational hybrid with predominantly Blue-winged Warbler ancestry that happens to be homozygous recessive for the black throat allele, (2) Brewster’s hybrids can have majority ancestry from either parental species and are not homozygous recessive for the black throat allele, and (3) the “intergrades” we sequenced are multigenerational hybrids with predominantly Golden-winged Warbler ancestry that are homozygous recessive for the black throat allele. In other words, any categorization of hybrids into discrete types is based primarily on one discrete trait (throat patch), whereas there is continuous variation in carotenoid-based pigmentation and considerable underlying genomic variation in hybrids, particularly for Brewster’s individuals.

We note that even with genomic data, it is difficult to trace the complex ancestral history for multigenerational hybrids (Anderson and Thompson 2002, Milne and Abbott 2008), which may result from backcrosses (e.g., F2 x Golden-winged Warbler) or crosses between other hybrids with complex ancestry (e.g., F1 x F3, or F2 x backcrossed Golden-winged Warbler). Due to this, and because we are limited by our small sample size in this study, we do not assign individuals to specific hybrid classes and limit our inference of individual admixture to majority ancestry (hybrid index) and heterozygosity.

Underparts Color May Not Distinguish F1 Hybrids

Nichols (1908) hypothesized that F1 Vermivora hybrids would be Brewster’s type with white underparts, whereas Parkes (1951) hypothesized that F1s would be Brewster’s type with yellow underparts and that Brewster’s with white underparts are the product of backcrossing with Golden-winged Warblers. Our analyses of genomic admixture indicate that none of the hybrid individuals sequenced in this study are likely to be F1s (i.e. no individuals falling to the top of the triangle plot; Figure 1A). Due to the phenotypic variation they display, however, we can still address these hypotheses. We sequenced 6 Brewster’s hybrids that showed a wide range in expression of carotenoid-based coloration (Figure 1A). There are Brewster’s hybrids with yellow on their underparts that have hybrid index scores both >0.5 and <0.5 and the 2 Brewster’s hybrids with the whitest underparts have hybrid index scores <0.5, indicating predominantly Blue-winged Warbler ancestry. This result is inconsistent with Parkes’ (1951) hypothesis that Brewster’s hybrids with white underparts result from only backcrossing with Golden-winged Warblers, as several Brewster’s individuals with high hybrid index scores (consistent with high Golden-winged Warbler ancestry) had moderate yellow scores (Figure 1A).

None of the hybrids included in this study had ancestry values that clearly indicate they were F1s; therefore, our data cannot directly address whether F1s are predominantly white or yellow. However, a previous low-resolution genotyping assay by Toews et al. (2016) found that the only individual that was heterozygous across all 6 divergent SNPs was identified as a Brewster’s warbler, which also had extensive yellow coloration. While this is consistent with F1s having yellow underparts (i.e. Parkes 1951), additional genotype and phenotype data of individuals with known pedigrees will be desirable to explicitly test this hypothesis. Thus, while it may be possible to infer a Vermivora hybrid’s phenotype based on its genomic ancestry across the plumage-associated genes, vice versa (i.e. predicting a bird’s genetic ancestry based on only plumage) is much more challenging.

Further complicating the genotype-phenotypic associations, we observed within-individual variation in the degree of carotenoid-based pigmentation for Brewster’s hybrids that were studied across multiple years (Figure 2). Across years, there was a tendency for their plumage to decrease in carotenoid pigmentation. One possible explanation is that this results from disruption of the carotenoid-processing pathway due to hybrid dysfunction. Carotenoid-based plumage is a character that functions as an ornament in signaling male quality in wood warblers (Freeman-Gallant et al. 2010, 2011, Jones et al. 2017) and thus may play a role in hybrid fitness. Previous studies have found little evidence for reduced hybrid fitness in this system (Vallender et al. 2007), although individual hybrids have not been followed long term (but see Carter 1944), and the relationship between male reproductive success and degree of carotenoid pigmentation has not been assessed. Future investigations of hybrid fitness that consider individual age and variation in carotenoid pigmentation will provide more insight on this matter.

Alternately, this carotenoid-based plumage change could be due to seasonal shifts, or a result of age- or diet-based differences across years. Seasonal variation in carotenoid-based pigmentation has been observed in House Finches (Haemorhous mexicanus) and Great Tits (Parus major) (McGraw and Hill 2004, Figuerola and Senar 2005), likely due to feather or pigment degradation. In House Finches, adults express more carotenoids than hatch-year birds (Hill 1992, Inouye et al. 2001). In Golden-winged Warblers, first winter birds often have some yellow on their underparts but adults have white underparts (Parkes 1951). This reduction in carotenoid pigmentation is consistent with what we observed here for Brewster’s hybrids, but in the current study we aged the hybrids to be “after second year” upon first capture. Thus, transition from hatch-year to adult plumage cannot explain our observations. We are aware of no other studies that have followed annual variation in Vermivora hybrid plumage patterns, aside from Carter (1944), who spent 6 yr observing the same Brewster’s warbler hybrid from 1922 to 1927. However, while Carter (1944) mentioned other aspects of annual variation in this bird’s feather patterns (e.g., a head scar that developed between 1924 and 1925), he did not note any variation in carotenoid-based plumage.

In captive feeding experiments, Hill (1992, 1993) showed that variation in plumage coloration for adult House Finches varies strongly as a function of dietary access to carotenoids during molt as coloration tends to converge for birds given the same diet. Further, carotenoid pigmentation was consistent across years for adults (Hill 1992). In wood warblers, environmental effects on pigmentation have not been studied, but our results may suggest an interaction effect of age and/or diet with genetics on carotenoid pigmentation in Vermivora warblers.

Taken together, our results suggest that underparts color is a complex trait with a genetic basis that is not straightforward and, thus, it may not be under the control of a single gene as previously predicted (Nichols 1908, Parkes 1951). Further, because our results indicate that Brewster’s hybrids beyond the first generation can display underparts ranging from white to yellow-washed and the degree of yellow can vary across years, we posit that it is not likely that F1 Brewster’s can be distinguished from multigenerational Brewster’s hybrids based on the degree of carotenoid pigmentation on their underparts.

Variation Near ASIP Promotor is Linked to Throat Patch Color

Our analyses of genomic variation and throat patch color confirms that this phenotype is likely under the control of a small region on warbler scaffold 299 (Zebra Finch chromosome 20). This is consistent with a previous study that discovered low levels of genome-wide divergence between Blue-winged Warblers and Golden-winged Warblers, but high differentiation near ASIP (Toews et al. 2016). The previous study did not include re-sequencing data from hybrids but, by doing so here, we narrowed the candidate throat patch region from >50 Kb down to ~10–15 Kb. First, on the candidate warbler scaffold, we observed greater nucleotide diversity around position 390 Kb for Brewster’s hybrids than for the other plumage types, which would be expected to be homozygous for the throat patch allele if it is controlled by a single locus inherited in Mendelian fashion. Second, in our iterative principle component analysis, we observed a switch in the position of the Lawrence’s hybrid (the only hybrid in our analysis with a black throat patch, but otherwise yellow Blue-winged Warbler-like plumage) from clustering intermediate to the parental species before position 385 Kb and after position 400 Kb, but clustering with Golden-winged Warblers and intergrades (which express the black throat) within this window.

The region between 385 and 400 Kb occurs in an area of several SNPs with nearly fixed differences between Blue-winged and Golden-winged Warblers upstream of the ASIP gene and, notably, upstream of its predicted promoter based on alignment of chicken sequences (Toews et al. 2016). This may suggest that wood warblers have different or additional promotor regions than chickens to regulate ASIP. Genetic differentiation in sequence near ASIP has also been linked to species differences in plumage pigmentation in Sporophila seedeaters (Campagna et al. 2017), in Lonchura munias (Stryjewski and Sorenson 2017), and in Setophaga warblers (Wang et al. 2019). Similarly, ASIP expression contributes to differences in feather pigmentation in juncos (Abolins-Abols et al. 2018). Additional studies are necessary to test whether divergence in the ~10–15 Kb candidate region we identified here underlies expression differences in ASIP, and whether throat patch color varies as a result of this in Vermivora and in other wood warbler species. It will be possible to test this hypothesis using gene expression analyses such as RNAseq (e.g., Abolins-Abols et al. 2018).

Although our analysis includes few hybrids and only one Lawrence’s individual, we were able to identify a small genomic region that likely underlies the discrete throat patch phenotype in this system. Expanded sequencing of Lawrence’s hybrids would be desirable to replicate our results, but this may be difficult since hybrids of this type are rare. Notably, the Lawrence’s hybrid sequenced here is not the same individual genotyped in the previous study linking ASIP to the throat patch phenotype (Toews et al. 2016). Thus, our results independently implicate this as a candidate region and highlight the power of sequencing natural hybrids to link discrete phenotypes to genomic variation.

CONCLUSIONS

Here, we used whole genome re-sequencing to explore links between genomic and phenotypic variation in hybrids of Blue-winged Warblers and Golden-winged Warblers. Our results indicate that hybrid types are quite variable and that F1 hybrids may not be easily distinguished by coloration of their underparts, as previously predicted (Nichols 1908, Parkes 1951). Further, it seems that underparts coloration does not have a straightforward genomic basis, unlike throat patch coloration, which we have linked to an ~10–15 Kb region upstream of the ASIP gene. Our observations also suggest there may be an interaction between genetics and diet or age on carotenoid-based pigmentation in hybrids, and future studies that include hybrids should consider that it is not necessarily a fixed trait.

APPENDIX TABLE 2. Sampling localities for individuals included in this study. BWWA = Blue-winged Warbler, GWWA = Golden-winged Warbler. State: NY = New York, MI = Michigan, PA = Pennsylvania, WI = Wisconsin. See Figure 1B in the main text for illustrations of plumage types.

Plumage type	ID	Locality	State	Latitude	Longitude
BWWA	215081703	Harriman St Park	NY	41.233434	–74.091252
BWWA	215081704	Sterling Forest	NY	41.234886	–74.237887
BWWA	215081705	Sterling Forest	NY	41.225508	–74.213411
BWWA	215081708	Sterling Forest	NY	41.188546	–74.260289
BWWA	215081709	Sterling Forest	NY	41.189318	–74.260195
BWWA	215081720	Rossie	NY	44.390593	–75.698829
BWWA	215081724	Rossie	NY	44.400149	–75.695038
BWWA	215081741	Finger Lakes	NY	42.482955	–76.790231
BWWA	215081742	Warren Road	NY	42.493499	–76.449233
BWWA	215081743	Summer Hill State Forest	NY	42.664890	–76.346280
GWWA	215081701	Sterling Forest	NY	41.210971	–74.236410
GWWA	215081707	Sterling Forest	NY	41.197464	–74.254768
GWWA	215081711	Southfields	NY	41.261847	–74.174682
GWWA	215081712	Southfields	NY	41.262409	–74.174583
GWWA	215081714	Rossie	NY	44.375757	–75.675113
GWWA	215081716	Rossie	NY	44.373562	–75.694552
GWWA	215081721	Rossie	NY	44.372909	–75.704105
GWWA	215081725	Lonesome Bay State Forest	NY	44.419696	–75.667912
GWWA	215081733	Rossie	NY	44.388354	–75.638707
GWWA	215081738	Rossie	NY	44.380308	–75.598252
Brewster’s	215081706	Sterling Forest	NY	41.196066	–74.255980
Brewster’s	215081713	Rossie	NY	44.383480	–75.661014
Brewster’s	MIH02	Secord	MI	44.013112	–84.307888
Brewster’s	PAH01	Bald Eagle State Park	PA	41.008718	–77.676564
Brewster’s	PAH06	Stone Valley Forest	PA	40.666372	–77.934267
Brewster’s	WIH03	Wood County	WI	44.309842	–90.177886
Intergrade	215081718	Rossie	NY	44.384973	–75.700134
Intergrade	215081729	Lonesome Bay State Forest	NY	44.426426	–75.663025
Lawrence’s	PE25T03	Long Meadow Drive	NY	41.183438	–74.256229

APPENDIX TABLE 3. Plumage scores for individuals sequenced in this study. BWWA = Blue-winged Warbler, GWWA = Golden-winged Warbler, Brewster’s = Brewster’s hybrid, Intergrade = backcrossed hybrid, Lawrence’s = Lawrence’s hybrid, wbar = wingbar, sum plumage = sum of scores for all plumage characters (gray columns), sum yellow = sum of scores for only carotenoid-based plumage characters (gray columns with asterisk in heading). We report the “modified sum yellow score” in the main text, which is the inverted “sum yellow” score, calculated as the sum yellow score subtracted from the maximum (i.e. 31-sum yellow). Higher modified sum yellow scores indicate plumage that is more yellow overall. See Figure 1B in the main text for illustrations of plumage types.

Plumage type	ID	wbar color*	wbar width	Throat	Nape*	Back*	Rump*	Breast*	Belly*	Mask	Eye- line*	Mustache*	Sum plumage (out of 38)	Sum yellow (out of 31)	Modified sum yellow score (31-sum yellow)
BWWA	215081703	2	0	0	0	0	0	0	0	0	0	0	2	2	29
BWWA	215081704	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081705	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081708	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081709	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081720	0	0	0	0	0	0	0	0	0	0	0	0	0	31
BWWA	215081724	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081741	1	0	0	0	0	0	0	0	0	0	0	1	1	30
BWWA	215081742	0	0	0	0	0	0	0	0	0	0	0	0	0	31
BWWA	215081743	1	0	0	0	0	0	0	0	0	0	0	1	1	30
GWWA	215081701	5	3	2	4	4	4	4	4	2	3	3	38	31	0
GWWA	215081707	5	3	2	4	4	4	4	4	2	1	3	36	29	2
GWWA	215081711	4	2	2	4	4	4	4	4	2	3	3	36	30	1
GWWA	215081712	5	3	2	4	4	4	4	4	2	3	3	38	31	0
GWWA	215081714	5	3	2	4	4	4	4	4	2	3	3	38	31	0
GWWA	215081716	5	2	2	4	4	4	4	4	2	2	3	36	30	1
GWWA	215081721	5	2	2	4	4	4	4	4	2	3	3	37	31	0
GWWA	215081725	4	3	2	4	4	4	4	4	2	2	3	36	29	2
GWWA	215081733	5	2	2	4	4	4	4	4	2	3	3	37	31	0
GWWA	215081738	4	2	2	4	4	4	4	4	2	3	3	36	30	1
Brewster’s	215081706	4	2	0	3	3	3	2	2	0	2	2	23	21	10
Brewster’s	215081713	4	1	0	3	2	2	3	4	0	1	2	22	21	10
Brewster’s a	MIH02_2015	4	1	0	1	2	2	0	1	0	1	2	14	13	18
Brewster’s a	MIH02_2016	4	1	0	3	3	3	2	3	0	1	3	23	22	9
Brewster’s a	PAH01_2015	5	3	0	1	2	2	2	1	0	3	3	22	19	12
Brewster’s a	PAH01_2016	4	2	0	3	4	4	3	4	0	2	3	29	27	4
Brewster’s	PAH06_2015	5	2	0	2	3	3	2	3	0	1	2	23	21	10
Brewster’s a	WIH03_2015	4	1	0	2	4	4	3	3	0	2	3	26	25	6
Brewster’s a	WIH03_2016	4	2	0	3	4	4	3	3	0	2	3	28	26	5
Brewster’s a	WIH03_2017	4	2	0	4	4	4	3	3	0	2	3	29	27	4
Intergrade	215081718	3	2	2	3	2	2	1	3	2	1	2	23	17	14
Intergrade	215081729	5	3	1	3	3	3	2	3	2	2	3	30	24	7
Lawrence’s	PE25T03	2	0	1	0	0	0	0	0	2	0	0	5	2	29

^a Indicates score for individuals studied in multiple years (year reported is appended to individual ID in column 2).

APPENDIX TABLE 4. Number of genomic sites included in our analyses. Genotype likelihood approach indicates we used ANGSD to estimate genotype likelihoods for analyses. Genotype call approach indicates we used HaplotypeCaller to call genotypes for analyses.

Warbler scaffold:position (bp)	N sites, genotype likelihood approach, principle component analyses	N sites, genotype call approach, hybrid admixture analyses	N sites, genotype call approach, nucleotide diversity analyses
120:400,000–450,000	4,270	91	NA
24:8,600,000–8,880,000	12,287	338	NA
Warbler scaffold 299 a	NA	NA	30,747
299:372,712–394,284	2,033	54	NA
299:380,000–385,000 b	542	NA	NA
299:385,000–390,000 b	394	NA	NA
299:390,000–395,000 b	561	NA	NA
299:395,000–400,000 b	426	NA	NA
299:400,000–405,000 b	627	NA	NA
299:405,000–410,000 b	555	NA	NA
38:2,400,000–2,870,000	23,077	413	NA
563:42,960–52,963	586	16	NA
653:92,441–105,000	845	48	NA

^a Indicates the entire scaffold was used in analysis, otherwise positions denoted were within divergent regions in Toews et al. (2016).

^b Denotes 5 Kb windows included in our iterative PCA analysis of warbler scaffold 299.

APPENDIX FIGURE 5.

Close relationship between admixture scores calculated in this study using 960 SNPs (Introgress score), and in a previous RFLP genotyping assay (Toews et al. 2016) using 6 SNPs (RFLP score) for the 4 hybrid individuals that were included in both studies. The gray line represents a 1:1 relationship.

ACKNOWLEDGMENTS

We thank J. D. Curlis and G. Hill for helpful discussion, B. G. Butcher for assistance in lab work, and N. Hofmeister for field sampling.

Funding statement: Field sampling and sequencing were supported by the Cornell Lab of Ornithology, U.S. Fish and Wildlife Service, U.S. Geological Survey, and the National Science Foundation (award #1202729 to H.M.S.).

Ethics statement: Work was conducted under permits from the U.S. Fish and Wildlife Service (permit #24034 to D.P.L.T.; permit #21631 to D. E. Andersen), the New York State Department of Environment and Conservation (permit #153 to D.P.L.T.), Michigan Department of Natural Resources, Pennsylvania Game Commission, and Wisconsin Department of Natural Resources. All animal research was conducted ethically under the purview of the Cornell University Institutional Animal Care and Use Committee (IACUC; protocol #2015-0065 to I.J.L.), and University of Minnesota IACUC (protocol #1302-30370A).

Author contributions: M.D.B. and D.P.L.T. conceived the study, designed the methods, and wrote the paper, G.R.K., H.M.S., S.A.T., and D.P.L.T. collected samples, M.D.B. analyzed the data, and I.J.L. and H.M.S. contributed substantial resources and funding. All authors contributed to the editing of the manuscript and approved the final version.

Data depository: Re-sequencing data for all parental individuals, Brewster’s and intergrade hybrids are available at the NCBI Short Read Archive under BioProject PRJNA325126. Digital photographs of all individuals sequenced in this study are available in the Dryad dataset at https://doi.org/10.5061/dryad.tx95x69tk.