Comparative genomics of two Empidonax flycatchers reveal candidate genes for bird song production

Abstract

Whole-genome-level comparisons of sister taxa that vary in phenotype against a background of high genomic similarity can be used to identify the genomic regions that might underlie their phenotypic differences. In wild birds, this exploratory approach has detected markers associated with plumage coloration, beak and wing morphology, and complex behavioral traits like migration. Here, we use genomic comparisons of two closely related suboscine flycatchers (Empidonax difficilis and E. occidentalis) and their hybrids to search for candidate genes underlying their variation in innate vocal signals. We sequenced the genomes of 20 flycatchers that sang one of two species-specific pure song types and 14 putative hybrid individuals with intermediate song types. In the resulting genomic comparisons, we found six areas of high differentiation that may be associated with variation in nonlearned songs. These narrow regions of genomic differentiation contain a total of 67 described genes, of which three have been previously associated with forms of language impairment and dyslexia in humans and 18 are known to be differentially expressed in the song nuclei regions of the avian brain compared with adjacent parts of the avian brain. This “natural experiment” therefore may help identify loci associated with song differences that merit further study across bird lineages with both learned and innate vocalizations.

Introduction

Much of our current knowledge about the genomic basis of vocal communication in vertebrates derives from research on speech disorders in humans (Graham & Fisher, 2013), the only species capable of developing complex languages (Fitch, 2000). However, a key precursor for the development of human spoken language is the capacity for vocal imitation (Lattenkamp & Vernes, 2018), an ability present in other mammals, including some cetaceans, pinnipeds, elephants, bats, as well as in some birds, primarily parrots, hummingbirds, and songbirds (Janik & Slater, 1997; Petkov & Jarvis, 2012). As nonhuman primates and other animals used as model organisms for medical research are considered to lack the capacity for imitative vocal learning (but see Fischer & Hammerschmidt, 2020; Martins & Boeckx, 2020), research in this area has focused instead on a songbird, the Zebra Finch (Taenopygia guttata, Scharff & Adam, 2013), resulting in substantial advances in the genomics of vocal learning. In many groups of birds, vocal signals are diverse in form, but they are innate rather than learned. The order Passeriformes includes two suborders, the song-learning oscines or songbirds (Passeriformes: Passeri) and the innate-song suboscines (Passeriformes: Tyrannida). In suboscine species where vocal development has been studied experimentally, normal production of songs was achieved even in the absence of a tutor, indicating these signals are innate (Kroodsma, 1984; Touchton et al., 2014, but see ten Cate, 2021 for a review on the evidence of vocal learning in two suboscine species in the family Cotingidae). Comparative genomics of such species could provide the means to investigate the underlying genetic architecture of nonlearned vocal communication.

In taxa that are not amenable to laboratory-based studies, comparisons among closely related species can provide a form of “natural experiment” that allows for the identification of candidates for the genetic basis of variable traits. Whole-genome comparisons of extremely similar taxa with low levels of genetic differentiation may help find the genomic regions that underlie their few phenotypic differences, which are expected to show divergence levels that stand out prominently against a low background of differentiation across the genome (Seehausen et al., 2014). In birds, this comparative approach has proven useful for detecting markers related to morphological phenotypes such as plumage coloration (Campagna et al., 2017; Knief et al., 2019; Toews et al., 2016), wing morphology (Campagna et al., 2019), beak morphology (Chaves et al., 2016; Lamichhaney et al., 2016; Lawson & Petren, 2017), and behavioral traits such as migration (Merlin & Liedvogel, 2019; Ruegg et al., 2014a, 2014b; Toews et al., 2019). Similarly, variation in innate vocal signals might be linked to specific genomic variation between closely related species. Song differences have been shown to accumulate relatively quickly between diverging species with innate songs, even faster than between vocal-learning birds (Freeman et al., 2017).

Kroodsma (1984) provided the first direct experimental evidence that song development in suboscine passerines (Passeriformes: Tyrannida) is innate by working with two tyrannid flycatcher species in the genus Empidonax, willow flycatchers (E. trailli), and alder flycatchers (E. alnorum). When individuals of these two species are removed from the nest and tutored with either heterospecific or conspecific songs, all individuals ultimately produced the normal song of their species regardless of the tutelage to which they were exposed (Kroodsma, 1984). Here, we investigate whole-genome differences between two sister species in the same genus, Pacific-slope and Cordilleran flycatchers (Empidonax difficilis and E. occidentalis, Figure 1A). These two taxa were long considered conspecific but were taxonomically split into different species (Monroe et al., 1989) due primarily to the discovery of differences in their innate vocalizations (as they are virtually indistinguishable based on plumage coloration, morphology, and other physical traits, Lowther et al., 2020). As song is an important reproductive isolating mechanism for birds (Uy et al., 2018), the evolutionary divergence of bird song is relevant to speciation, either as a product of it or as a driving force. In one scenario, song divergence is a byproduct of genetic differentiation accumulated over prolonged periods of allopatry, stabilizing via reinforcement into distinct song types when populations later come into secondary contact (Butlin, 1987; Dobzhansky, 1937; Servedio & Noor, 2003). In another scenario, signals diverge early in differentiation (driven by sexual selection and/or ecological adaptation and/or drift), and the resulting premating isolation based on assortative mating by song precedes genome-wide differentiation (Coyne, 1992; West-Eberhard, 1983). Therefore, when comparing the genomes of two recently diverged flycatcher species with diagnostic song differences, we expect to find areas of the genome containing either loci that are associated with gene flow resistance or loci that have been subject to recent selection.

Empidonax difficilis and E. occidentalis have a distinct but shallow mtDNA divergence and likely diverged approximately 350,000 years ago (0.7% mitochondrial distance; Johnson & Cicero, 2002). Initially, it was believed that their only area of sympatry was in the Siskiyou Mountains in northeastern California where the two taxa function as reproductively isolated species (Johnson, 1980, 1994; Johnson & Marten, 1988). However, a more recent study incorporating populations from interior southwestern Canada found evidence of ongoing hybridization (Rush et al., 2009), consistent with reports of flycatchers in this region having vocal features that are intermediate between the standard Pacific-slope and Cordilleran types. Noteworthy, northern E. occidentalis populations (from the United States) are more closely related to E. difficilis than to E. occidentalis populations from Mexico (Linck et al., 2019). Therefore, E. occidentalis as currently defined may be paraphyletic with respect to E. difficilis, as individuals of both species from the United States grouped in sister clades nested within E. occidentalis from Mexico (Linck et al., 2019). The diagnosable innate song differences between these species, their shallow overall divergence, and the availability of hybrids with intermediate song types make this an intriguing system to prospect for candidate areas of the genome related to innate vocal signals.

Methods

Sampling

Pacific-slope and Cordilleran flycatchers are distributed across much of temperate western North America and Mexico, with several known areas of secondary contact (Johnson, 1980; Linck et al., 2019; Rush et al., 2009). We selected tissue samples from the collection of the Museum of Vertebrate Zoology at the University of California, Berkeley, corresponding to 34 individuals that were collected from 20 localities across the northern populations of these species (Figure 2A), avoiding E. difficilis individuals from Mexico (Linck et al., 2019). The individuals were recorded prior to collection; therefore, each tissue sample is directly linked to song recordings from the same bird (Supplementary Table S1). When these individuals were collected, those identified as potential migrants (displaying high levels of subcutaneous fat and the absence of enlarged gonads) were noted and excluded from selection in this study (Rush, 2014). We classified individuals into three song groups (Pacific-slope, Cordilleran, or intermediate, see below) based on a previously obtained “song score” that distinguishes pure and intermediate song types and was based on the songs of 608 individuals from 68 sites (see Table A4 in Rush, 2014), including those sequenced here. Briefly, the songs of these flycatchers are comprised of three distinct elements that are most often delivered in a repetitive 1–2–3 order (Figure 1A). The score was obtained by running a principal components analysis (PCA) focused on the second element of the song, since it is the most structurally complex and best represents differences between Pacific-slope and Cordilleran song types. The acoustic variables measured were peak frequency, duration, sharpness of the frequency peak (i.e., the absolute value of the slope from the frequency at peak frequency minus 10 ms to peak frequency + the slope from the peak frequency to the frequency at peak frequency plus 10 ms), the proportion of the entire song comprised of the second part of the song (i.e., after the lowest inflection point), the change in frequency from the peak to the lowest inflection point, and the presence or absence of an amplitude gap at the lowest inflection point. For further details, see Rush (2014). The first PC (Song PC1) axis of song variation, with an eigenvalue of 3.56, explained 51% of the variation and separated individuals from the parental populations into two largely nonoverlapping clusters (see Chapter 2, figure 7 in Rush, 2014). PC1 values were used to classify the flycatchers sequenced for this study into three groups (Figure 1B): those with a song score > 1 were assigned to a Pacific-slope song group (n = 10), those with a song score < −1 were assigned to the Cordilleran song group (n = 10), and those with intermediate values were considered to have an intermediate song (n = 14) (Supplementary Table S1). We refer hereafter to these groups as “PSFL-SONG,” “COFL-SONG,” and “INTER-SONG,” respectively.

Sequencing, alignment to reference genome, and SNP discovery

We used the Puregene tissue extraction protocol (Gentra Systems) to purify DNA samples and prepared individually barcoded libraries following the TruSeq Nano DNA library preparation kit protocol for an insert size of 350 bp (median insert size obtained: 314–442 bp). The 34 libraries were sequenced on three Illumina NextSeq 500 lanes at the Cornell Institute for Biotechnology resulting in an estimated average coverage ± standard deviation per individual of 6.2× ± 1.7.

We assessed the quality of individual libraries using fastqc version 0.11.5 (www.bioinformatics.babraham.ac.uk/projects/fastqc) and then removed adapters and performed quality filtering with AdapterRemoval 2.1.1 (Schubert et al., 2016). We allowed a minimum Phred quality score of 10 and merged overlapping paired-end reads. We used Chromosemble to align the willow flycatcher (Empidonax trailli) reference genome (GenBank assembly accession: GCA_003031625.1, Ruegg et al., 2018) to the Zebra Finch genome (GenBank assembly accession: GCA_003957565.2) to orient scaffolds and assemble them into chromosomes. We subsequently mapped the filtered reads to the improved E. trailli genome with the very sensitive local option implemented in Bowtie 2 2.3.4.3 (Langmead & Salzberg, 2012). The average mapping rate ± standard deviation across all samples was 97.4% ± 0.9. We converted sam files into bam format and subsequently sorted and indexed these files with SAMtools 1.9 (Li et al., 2009). We marked PCR duplicates with Picard Tools version 2.8.2 (https://broadinstitute.github.io/picard/) and performed SNP variant discovery and genotyping with the haplotype caller module in GATK 3.8.1 (McKenna et al., 2010). We compiled two vcf files, one comprising all the genotyped individuals (34 samples) and one with the individuals from each of the pure song type groups, PSFL-SONG and COFL-SONG (20 samples). For each vcf file, we removed variants based on the following hard filtering parameters: QD < 2, FS > 40.0, MQ < 20.0, and HaplotypeScore > 12.0 using GATK 3.8.1. We then used VCFtools v. 0.01.16 (Danecek et al., 2011) to filter out nonbiallelic variants, with a minor allele count smaller than 3, with mean depth of coverage smaller than 3 or greater than 50, and/or with >20% missing data across the data set. This pipeline produced 13,587,433 SNPs across all individuals (average depth of coverage 5.9× ± 1.6×) and 10,552,562 SNPs in the subset of pure song individuals (average depth of coverage 6.2× ± 1.9×).

Population genomics analyses

We first ordinated the SNP variation across all 34 samples in a PCA via the package SNPRelate version 3.3 (Zheng et al., 2012) in R version 3.6.1 (R Core Team, 2019). We also ran an ancestry-estimation analysis in Admixture 1.3 (Alexander et al., 2009) at K = 2, to estimate the levels of ancestry of individuals singing pure and intermediate song types. To avoid including SNPs in high linkage disequilibrium we first thinned the vcf file, keeping sites separated by a minimum of 1,000 base pairs.

We then used VCFtools to search for divergent areas of the genome by calculating F_ST values between PSFL-SONG and COFL-SONG for nonoverlapping 25-kbp windows with at least 10 SNPs and for individual SNPs. VCFtools calculates the F_ST estimator proposed by Weir and Cockerham (1984) and also provides weighted averages for the calculations over windows of loci, as recommended by Weir and Cockerham (1984). We considered outliers any two or more consecutive windows with a weighted F_ST higher than 0.5 and containing at least one SNP near fixation (F_ST > 0.85). These criteria focused on the main divergence peaks between the species. None of the scaffolds that failed to align with the Zebra Finch chromosomes contained outlier windows. To search for possible genes of interest, we delimited areas containing the outlier windows plus 50 kbp to either side (Supplementary Table S2) and used BLAST (Altschul et al., 1990) to align them to the E. trailli genome (e-value < 10⁻¹⁰). This also allowed us to confirm that the outlier areas we found (using the E. trailli genome realigned to the Zebra Finch genome as the reference) corresponded to only one or two scaffolds from the original (unordered) E. trailli reference genome.

We then focused on areas containing the outlier windows plus 1,000 kbp on either side (Supplementary Table S2) and calculated weighted F_ST values for nonoverlapping 5-kbp windows with at least 2 SNPs, along with D_XY, Tajima’s D, and nucleotide diversity (π). We calculated Tajima’s D using only variant sites (i.e., SNPs) as for the F_ST analyses, but we first removed the minimum allele count filter to avoid biases toward positive values that result from excluding alleles segregating at lower frequencies. We used all sequenced sites, both variant and invariant, to calculate D_XY and π. To obtain these, we also used the GATK command GenotypeGVCFs to compile vcf files containing the sequences for the 20 PSFL-SONG and COFL-SONG individuals, this time using the option “-allSites.” As vcf files including all sites (not just SNPs) can be extremely large, we compiled one vcf per chromosome. To remove sites with MQ < 20 and avoid poorly mapped reads, we used a custom-written script by Irwin et al. (2016) to filter out sites with MQ lower than 20.0 but leave in invariant sites without an MQ score (which would have been removed had we been using GATK as we did previously). We used VCFtools to filter sites with mean depth of coverage smaller than 3 or greater than 50, and/or with >20% missing data across the data set. D_XY and π values were calculated using custom-written scripts by S. H. Martin (downloaded from https://github.com/simonhmartin/genomics_general/). We used the package qqman 0.1.4 (Turner, 2017) in R to build a Manhattan plot.

We also used the SNPRelate package (Zheng et al., 2012) to run principal component analyses for each chromosome using the nearly fixed SNPs of the outlier windows (F_ST > 0.85). We used plink version 1.9 (Purcell et al., 2007) to calculate the r² statistic in the outlier areas to explore patterns of linkage disequilibrium. We considered PSFL-SONG, COFL-SONG individuals and INTER-SONG individuals separately. We first estimated r² using the nearly fixed SNPs across all of the outlier windows together, and we then did the same for each outlier area in each chromosome (±100 kbp, see Supplementary Table S2), considering SNPs with a minimum allele frequency of 0.25.

Additionally, we analyzed our data including both pure and intermediate song individuals in the context of a genome-wide association analysis in Gemma version 0.98.4 (Zhou & Stephens, 2014). This analysis attempts to fit univariate linear mixed models using the Song PC1 as a phenotype and including an interindividual relatedness matrix among all samples as a covariate to account for potential population structure. However, given our sample size, we were underpowered to detect associations between genotypes and phenotypes using this strategy, as we did not find regions showing statistical association values above the threshold of significance established after Bonferroni correcting for multiple comparisons (~13 million SNPs). We then used VCFtools to reduce the number of SNPs, and therefore the number of comparisons, by keeping only those with no missing information (~3 million SNPs); we still found no regions with statistical association levels above the Bonferroni-corrected threshold (Supplementary Figure S1).

Results

We used a PCA to summarize the genetic relationships between 34 flycatcher individuals based on ~13.6 million genome-wide SNPs (Figure 2B) and found that the first and second PCA axes collectively capture almost 10% of the total genomic variance. PSFL-SONG and COFL-SONG individuals with typical songs do not overlap in the genetic PCA, although some individuals fall close to individuals classified via song to the other species (Figure 2B). INTER-SONG individuals (likely hybrids with intermediate songs) primarily overlap genomically with both pure song groups, a pattern consistent with the presence of extensive backcrossing into either genomic background. The degree of overlap (and also perhaps of backcrossing) is higher with PSFL-SONG, which is also true geographically in our sample of INTER-SONG individuals. Consistently, our ancestry-estimation analyses (Figure 2C) showed that while all but one COFL-SONG individuals are not admixed, a few PSFL-SONG individuals have partial ancestry from the other species (Figure 2C). Most INTER-SONG individuals have varying levels of ancestry, suggesting the INTER-SONG group is primarily made up of a mix of F1 and backcrossed individuals. Genomic PC1 and Song PC1 show a significant, positive correlation (R² = 0.8, P < 0.001; Figure 2D).

The Manhattan plot of whole-genome differentiation between PSFL-SONG and COFL-SONG (pure song types) illustrates their generally low background differentiation (weighted average F_ST = 0.067), with higher differentiation in the Z chromosome (weighted average F_ST = 0.132), consistent with previous results showing genetic differences accumulate faster in this chromosome compared to the rest of the genome (Irwin, 2018). A few notable peaks of elevated differentiation are also evident in this plot (Figure 3A). Across the genome, we identified 71 outlier windows, which altogether contain 452 of the total 1,273 SNPs near fixation between PSFL-SONG and COFL-SONG (F_ST > 0.85, Figure 3A). We also ran a genome-wide association analysis but did not have sufficient statistical power to find SNPs with association values above the threshold of significance established after applying a Bonferroni correction (−log10(p) = 7.8 when comparing across ~3 million SNPs with no missing information; Supplementary Figure S1). A total of 29 SNPs showed a p-value of association above an arbitrary threshold of −log10(p) = 5, of which 10 fell within the previously established outlier windows on chromosomes 1A and 10.

To better understand what processes could have generated these areas of elevated differentiation, we calculated D_XY, π, and Tajima’s D and explored the patterns of these statistics in the differentiated outlier areas between PSFL-SONG and COFL-SONG and the areas immediately surrounding them. D_XY is a measure of the absolute genetic differentiation that has accumulated over time between two populations, while F_ST is a measure of relative genetic differentiation. While areas of elevated F_ST coincide only partially with areas of elevated D_XY and lower π compared to the surrounding areas on chromosomes 1A and Z (Figure 3B and C), D_XY and F_ST are not significantly correlated in these areas for any chromosome (p > .05 for all cases, Supplementary Figure S2a). Nucleotide diversity (π) is very similar in both groups (Figure 3C). We averaged π values for PSFL-SONG and COFL-SONG and found that D_XY and average π are almost perfectly correlated across each chromosome (p < 2.2 × 10⁻¹⁶ and adjusted R² > 0.98 for all cases, Supplementary Figure S2b). This finding suggests that D_XY values are similar to the levels of genetic diversity in each species and that the genetic differentiation accumulated between the taxa is small, as expected at early stages of divergence (Burri, 2017). Calculations of Tajima’s D recover a variety of patterns across the F_ST outlier regions (Figure 3D). In most cases, both PSFL-SONG and COFL-SONG exhibit negative Tajima’s D values in the outlier region, which can indicate a recent selective sweep. Along some portions of chromosomes 1A, 2 and Z, values of Tajima’s D are slightly positive for COFL-SONG but have negative values within the outlier region for PSFL-SONG, suggesting that there could have been selection within the outlier region in the latter lineage.

We then aligned the outlier areas ± 50 kbp against the annotated E. trailli reference genome and found 67 genes (see Supplementary Tables S3 and S4) and 22 uncharacterized loci. Notably, three of the genes we recovered have associations with forms of speech impairment in humans: ERC1 in the chromosome 1A peak (apraxia of speech in children, Thevenon et al., 2013) and KIAA0319 and DCDC2 in the chromosome 2 peak (dyslexia, Scerri & Schulte-Körne, 2010). Some of the nearly fixed SNPs fall immediately downstream of the DCDC2 gene or within the KIAA0319 gene (Figure 4A), including the SNP with highest F_ST in this region, for which we show the genotypes for the sequenced individuals (Figure 4B). In songbirds, neurons dedicated to song learning and development are organized in particular clusters or nuclei within the avian brain (Nottebohm, 2005), and our gene list includes 18 loci that are differentially expressed in these nuclei compared to adjacent areas of the brain (CCDC77, CD247, DCDC2, ERC1, MAP1A, MPC2, MPZL1, MRS2, NINJ2, POU2F1, PTPRD, RASGRF1, TUBGCP4, WNK1, KANK1, SLC16A3, FASN, and TMTC, see Supplementary Tables S3–S6 in Lovell et al., 2018). These genes are distributed mainly across the outlier areas on chromosomes 1, 1A, 2, and 10 (Supplementary Table S3). Vocal differences could also be a consequence of morphological differences in the avian vocal organ (the syrinx, Garcia et al., 2017), or parts of the tract such as the trachea or the beak (Fitch, 1999; García & Tubaro, 2018; Podos, 2001). However, we did not find overlap between the genes in outlier regions and a list of genes that are differentially expressed in facial tissues from three different avian species (Brugmann et al., 2010) or with a list of candidate genes for beak length identified with a whole-genome resequencing approach in a goose species (Huang et al., 2022).

We ran further PC analyses on sites with F_ST > 0.85 for the outlier windows in each chromosome and calculated correlations between the Genomic PC1 and Song PC1 (Figure 5A). Genomic PC1 and Song PC1 also showed a significant, positive correlation (R² > 0.6, P < 0.001) in all cases. The INTER-SONG individuals fall mostly between COFL-SONG and PSFL-SONG, with some individuals falling close to or even overlapping with either pure song type, which also supports the hypothesis that individuals singing intermediate songs are mostly backcrossed hybrids. Finally, we explored the patterns of linkage disequilibrium for the nearly fixed SNPs across all outlier windows combined (F_ST > 0.85) in PSFL-SONG/COFL-SONG individuals and INTER-SONG individuals. We found that LD tends to be lower both within and among peaks in individuals with an intermediate song type, as expected for genetically admixed individuals (Figure 5B, Supplementary Figure S2c).

Discussion

Whole-genome comparisons of recently diverged sister taxa can be used as exploratory tools to identify the genomic regions that likely underlie their phenotypic differences. This approach has proven useful in discovering markers related to complex behavioral traits, such as migration, in non-model bird species that are not amenable to controlled cross-breeding experiments in captivity (Merlin & Liedvogel, 2019; Ruegg et al., 2014a, 2014b; Toews et al., 2019). In an analogous manner, variation in innate vocal signals must be associated with currently unknown genomic differences. Here, we compared the genomes of two recently diverged flycatcher species whose most prominent difference is in their innate vocalizations and where individuals that sing intermediate songs are genetically admixed. We find that these birds have an overall low level of genomic differentiation, except for a few areas that contain, among others, genes expressed in brain tissues and differentially expressed in song nuclei within the avian brain. These loci may be potential candidates for the genetic basis underlying innate song differentiation, but our approach does not allow us to discard associations to nonvocal phenotypic differences that are not detectable to us. We therefore consider that our strategy identified a list of candidate genes that mediate phenotypic differences, of which a subset is likely to be involved in song divergence.

Assessments of differentiation between populations across millions of SNPs may reveal significantly heterogeneous patterns of variation across the genome. In the initial stages of divergence, genetic differentiation is typically low across the genome except for a few areas; as time since divergence increases, these areas of elevated differentiation are expected to increase in both number and size (Flaxman et al., 2014). We compared the genomes of PSFL-SONG and COFL-SONG and found results consistent with recent divergence, as previously estimated based on mitochondrial DNA divergence at approximately 350,000 years (Johnson & Cicero, 2002; Linck et al., 2019). Multiple evolutionary processes could have caused these few areas of elevated relative differentiation (F_ST) in the PSFL-SONG versus COFL-SONG genomes. The elevated regions on chromosomes 1A and Z show partially elevated D_XY values and therefore fit the pattern expected for a “speciation island,” an area of the genome that is resistant to gene flow between two populations that could contain loci relevant to the speciation process (Burri, 2017). However, the areas of high relative differentiation on chromosomes 1, 2, 10, and 18 do not show elevated absolute differentiation, which may be consistent with a scenario of positive selection in allopatry or perhaps background selection (Irwin et al., 2018). It has further been proposed that areas of elevated F_ST, but not D_XY, could be generated through background selection eliminating deleterious mutations (Burri, 2017; Cruickshank & Hahn, 2014). However, a recent simulation study concluded that although background selection can elevate F_ST, the conditions under which this happens require unrealistically high mutation rates (Charlesworth et al., 1997; Matthey-Doret & Whitlock, 2019). F_ST values are most likely not significantly affected by locus-to-locus variation in the intensity of background selection if realistic parameters are used (Matthey-Doret & Whitlock, 2019). Additionally, at early stages of differentiation, D_XY and π are highly correlated (Riesch et al., 2017), as is the case for the scaffolds where we found F_ST peaks. At such early stages, D_XY values are related more to the levels of genetic diversity present at the beginning of the divergence process than to absolute genetic differentiation accumulated between the taxa (Burri, 2017; Riesch et al., 2017). Therefore, we are confident that the areas of elevated F_ST we found are not “false positives” resulting from background selection, but rather represent areas that either contain loci relevant to reducing gene flow or have been subject to recent selective sweeps in at least one species (Hejase et al., 2020).

Of the 67 loci we identified in outlier areas, 18 are among those that are differentially expressed in the song nuclei compared with other parts of the avian brain in the Zebra Finch (Lovell et al., 2018), a vocal learner. Neurons dedicated to song learning and production in songbirds are organized in clusters or nuclei that are roughly analogous to different layers of the mammalian cortex (Nottebohm, 2005; Reiner et al., 2004). These discrete brain nuclei and their connections are part of two main branches, the anterior and posterior pathways, due to their location within the bird telencephalon (Nottebohm, 2005). Suboscine birds like Empidonax sp. have been much less studied but are thought to apparently lack this neural organization (Gahr, 2000). However, two species of suboscine, the eastern phoebe (Sayornis phoebe) and the common scale-backed antbird (Willisornis poecilinotus), are known to have brain regions that are morphologically similar to one such nucleus (the RA center) of vocal-learning species, along with similar patterns of gene expression (De Lima et al., 2015; Liu et al., 2013). Therefore, our results are consistent with the expectation that there may be some shared molecular pathways between vocal learners and nonlearners.

We note that additional genes could be involved in the song differences between our focal species, yet we decided to focus on the peaks showing the strongest signals of differentiation in our data set. Vocal differences can arise due to multiple factors acting at different levels beyond neural control. For example, it has been shown in suboscines, which are mostly innate singers, that the syrinx has three sound sources (the vibrating membranes), as opposed to oscine passerines that learn to sing and have two sources of sound. This level of morphological specialization in innate singers may represent an alternative path to vocal learning in the evolution of avian acoustic diversity (Garcia et al., 2017). Similarly, the avian beak has been described as a “magic trait,” meaning it is a trait normally under ecological selection, which can also affect the production of a key sexual signal such as song, and therefore indirectly lead to nonrandom mating (Podos et al., 2013; Servedio et al., 2011). As a preliminary approach, we also searched for candidate genes involved in variation in beak morphology, but we did not find genes in the outlier areas that mediate this aspect of phenotype (we note that beak morphology is similar between our study species). This result does not eliminate the possibility that genes acting at levels other than the neuronal control of song could be involved in vocal differences between our study species.

Traits in distantly related species that are analogous or convergent at organismal and phenotypic levels can be generated by one or more genes or genetic networks that are homologous (“deep homology,” see Shubin et al., 1997). A clear example is FoxP2, the first gene to be directly associated with speech impairments in humans and thus involved in vocal production. Extensive research in a bird model species, the Zebra Finch, has shown that FoxP2 also plays an important role in avian song development (Heston & White, 2015; Scharff & Adam, 2013), even though vocal learning arose independently in humans and songbirds. Here, we found that three genes that have been related to language impairments and dyslexia in humans are among the loci potentially underlying the vocal differences between these recently diverged Empidonax species. While human dyslexia manifests primarily as difficulty in learning to read and spell, it is considered a language-related disorder due to underlying deficits in relevant aspects of cognitive processing, such as the manipulation of phonemes (Graham & Fisher, 2013; Ramus, 2013). Further studies are required to confirm a direct link between these loci and song variation in birds, which would extend the “deep homology” in human and avian vocalizations to other genes beyond FoxP2 and to nonvocal learner bird species.

In the case of the extensively studied FoxP2 gene, a single mutation is directly related to language impairments (i.e., monogenic model, Mountford & Newbury, 2018). However, in the genes, we recovered as outliers (KIAA00319, DCDC2, and ERC1), several mutations may confer susceptibility to some type of language impairment or dyslexia (Mountford & Newbury, 2018). KIAA00319 and DCDC2 are two genes most studied in relation to dyslexia, and both are involved in neuronal migration (Centanni, 2020). ERC1 encodes a protein that belongs to the family of RIM-binding proteins, which regulate neurotransmitter release (Wang et al., 2002). Nondeleterious mutations in these loci could contribute to shaping vocal differences in Empidonax flycatchers and other species. KIAA00319 and DCDC2 were shown to have comparable evolutionary histories in vocal-learning and nonlearning species in both mammals and birds (Mozzi et al., 2016), but there is no further information on the role of these and the other loci we recovered in vocal communication in nonhuman species. Additionally, the list of genes we found in the outlier regions of the genome includes 18 loci of potential interest in the regulation of vocal communication because they are differentially expressed in the brain song nuclei of the Zebra Finch compared with adjacent areas of the brain (Lovell et al., 2018). Further research on these loci could provide a much better understanding of the genetic architecture of vocal communication, which likely involves several genes and genetic pathways.

In summary, our results provide an exploratory window toward discovering the loci underlying differences in a behavioral trait that has been known for decades to have a genetic basis. Although our comparative studies of differentiated species have substantial inherent limitations, they take advantage of a powerful “natural experiment” underlaid by hundreds of thousands of years of differentiation, which cannot be replicated in any experimental setting. Our results identify loci that are potentially associated with song differentiation in non-model species with innate vocalizations and that merit broader investigation in this context. Our understanding of the genetic basis of vocal communication and its variation across individuals and species requires further studies looking into differential expression of candidate genes and into the effects of knocking down such genes in model animals, as has been done for FoxP2.

Data availability

Genomic data have been archived in GenBank (BioProject ID PRJNA780874). Scripts used to run analyses and create figures are available from https://github.com/ncg37/empidonax.

Author contributions

N.C.G., L.C., R.C.K.B., and I.J.L. conceived the idea and designed the study. A.C.R. collected tissue samples and song recordings. R.C.K.B. provided access to tissue samples and data. N.C.G. performed lab work. N.C.G. and L.C. performed bioinformatic analyses. N.C.G. wrote the manuscript with input from all authors.

Conflict of Interest: The authors declare no conflict of interest.

Acknowledgments

We thank Carla Cicero for her help providing access to the tissue samples and Bronwyn Butcher for her assistance with DNA extraction and library preparation. We also thank members of the Lovette Lab at the Cornell Lab of Ornithology for their useful comments and suggestions during the development of this project. This work was supported by the Fuller Evolutionary Biology Program at the Cornell Lab of Ornithology.

References