Gene flow during glacial habitat shifts facilitates character displacement in a Neotropical flycatcher radiation
© The Author(s). 2017
Received: 7 April 2017
Accepted: 10 August 2017
Published: 1 September 2017
Pleistocene climatic fluctuations are known to be an engine of biotic diversification at higher latitudes, but their impact on highly diverse tropical areas such as the Andes remains less well-documented. Specifically, while periods of global cooling may have led to fragmentation and differentiation at colder latitudes, they may – at the same time – have led to connectivity among insular patches of montane tropical habitat with unknown consequences on diversification. In the present study we utilized ~5.5 kb of DNA sequence data from eight nuclear loci and one mitochondrial gene alongside diagnostic morphological and bioacoustic markers to test the effects of Pleistocene climatic fluctuations on diversification in a complex of Andean tyrant-flycatchers of the genus Elaenia.
Population genetic and phylogenetic approaches coupled with coalescent simulations demonstrated disparate levels of gene flow between the taxon chilensis and two parapatric Elaenia taxa predominantly during the last glacial period but not thereafter, possibly on account of downward shifts of montane forest habitat linking the populations of adjacent ridges. Additionally, morphological and bioacoustic analyses revealed a distinct pattern of character displacement in coloration and vocal traits between the two sympatric taxa albiceps and pallatangae, which were characterized by a lack of gene flow.
Our study demonstrates that global periods of cooling are likely to have facilitated gene flow among Andean montane Elaenia flycatchers that are more isolated from one another during warm interglacial periods such as the present era. We also identify a hitherto overlooked case of plumage and vocal character displacement, underpinning the complexities of gene flow patterns caused by Pleistocene climate change across the Andes.
Evolutionary processes of diversification rely on environmental change to bring about shifts in the distributions of habitats and species. One of the greatest planetary engines of environmental change driving biotic diversification over the last 2,000,000 years has been the cyclical succession of periods of global cooling, the so called ‘ice ages’ [1–3]. Over the Quaternary period, ice ages have regularly led to significant earth-historic change through fluctuations in temperature, sea level and extensive glaciations.
The effects of ice ages on the biota of the northern hemispheric temperate regions are relatively well understood [4–8]. But, despite the wealth of research that has gone into the investigation of the effects of Pleistocene glaciations on biotic differentiation in northern temperate latitudes, much less is known about how periods of global cooling have impacted the most species-rich tropical regions of the world (for exceptions see, e.g., [9–11]).
One of the most influential theories on speciation, Jürgen Haffer’s ‘refuge theory’ , which was widely followed in the 1980s and 1990s, postulates that elevated levels of species diversity in tropical rainforests are a product of rainforest contraction and savannah expansion during planetary cycles. However, earth-historic evidence for lowland rainforest contractions during glacial periods has failed to materialize (e.g. ), and tropical zoologists are amassing ever more evidence in support of species ages in most tropical lowland rainforest vertebrates that would pre-date Quaternary global glaciations [13–18]. Meanwhile, the role of ice ages in the diversification of tropical montane fauna has received relatively less attention, with most recent speciation events mainly attributed to tectonic processes such as mountain build-up (e.g. [19–23]).
The Andes are arguably the most species-rich area on Earth, with levels of taxon diversity along their eastern rim considerably exceeding most other regions in the world for many animal classes [24, 25]. While the southern Andes in Chile, Peru, Argentina and Bolivia are high enough for limited glaciation, a large part of the Andes would have been unaffected by the advance of ice sheets during periods of global cooling. However, global ice ages may have substantially affected the biota of the Andes through global drops in temperature by up to 8 °C , which would have led to a downslope shift of various assemblages of montane forest habitat possibly by several hundred meters [27, 28]. Many montane species that are nowadays restricted to fragmented pockets of mountaintop habitat isolated from neighboring populations by valleys may have enjoyed increased connectivity during glacial periods when their montane forest would have linked up across lower-lying ridges. There is limited research that has addressed the role of Pleistocene climate change in potentially counteracting speciation by increasing connectivity among isolated populations, especially in such hyper-diverse areas as the Andes that are known to be dominated by species of young age [21, 29].
In this study, we investigate patterns of gene flow in a young Andean-centered radiation of tyrant-flycatchers, the so-called ‘montane Elaenia clade’ . In particular, we are interested in one species assemblage within the montane Elaenia clade, the so-called ‘albiceps complex’, that consists of three forms seemingly closely related to one another  but with differing levels of phenotypic and genetic divergence. Using a variety of population genetic and phylogenetic approaches, including coalescent simulations, and ~5.5 kb of DNA sequence data from dozens of samples from across the entire range of these species, we test whether recent gene flow among members of the ‘albiceps complex’, if any, is uniform across time or occurred disproportionately during the most recent period of global cooling. We also analyze whether patterns of gene flow support the hypothesis of character displacement having led to pronounced differences in phenotypic characters between two sympatric members of this complex, compared to less phenotypic differentiation in comparisons involving non-sympatric forms. Our study is one of few rigorous tests of the impacts of Pleistocene climate change on Neotropical montane biota that is based on a combination of acoustic, morphological and genetic data.
We compiled a total of 50 vocal recordings of albiceps, chilensis and pallatangae from public sound libraries, including xeno-canto (http://www.xeno-canto.org) and Macaulay Library (http://macaulaylibrary.org). Details of all recordings analyzed, including localities, dates, and names of recordists, are provided in Additional file 1: Table S1. All our 50 recordings are homologous with one another, constituting a monosyllabic call that is the most commonly-heard vocalization given by montane Elaenia species, not the more complicated polysyllabic ‘dawn song’ . This monosyllabic call note is distinctly different among the three Elaenia forms compared, but even so, we detected a limited number of obvious misidentifications that had been lodged with the sound libraries. We excluded recordings that may possibly be misidentified and ensured that the geographic locality of the recordings included in the analyses is in agreement with the newly assigned taxon’s distribution in every case. Using RAVEN PRO Version 1.5 (Bioacoustics Research Program, Cornell Laboratory of Ornithology, Ithaca, NY, USA), we measured seven vocal parameters: (i) call duration, (ii) proportion of time to reach maximum frequency, (iii) minimum frequency, (iv) maximum frequency, (v) peak frequency (defined as the frequency with the highest amplitude), (vi) ratio of bandwidth at maximum frequency versus at the start of the call, and (vii) ratio of bandwidth at maximum frequency versus at the end of the call (Additional file 1: Table S2). We used default settings of RAVEN PRO, except for window size, which was adjusted to 1024 to show optimal resolution of the spectrograms across all three taxa. We performed principal component analysis (PCA) on the vocal dataset to assess bioacoustic differences across the three taxa in R 3.2.1 .
To further highlight the difference in vocal traits, we calculated pairwise Euclidean distances between samples of different pairs of groups (albiceps-chilensis, chilensis-pallatangae and albiceps-pallatangae) by using their PCA coordinates from principal component 1 (PC1) and principal component 2 (PC2), as the latter two accounted for most of the variation within the dataset. We then used R to perform a Tukey’s test (as distances were normally distributed) to assess if these pairwise between-group distances are significantly different from one another.
DNA sampling, extraction and sequencing
In this study we sequenced nine loci (eight nuclear and one mitochondrial). Five of these (q4, q6, q8, q25, q26) are from among a panel of 26 anonymous loci distributed throughout the genome, primers for which were randomly generated using established protocols  and trialed for this project (Additional file 1: Table S4). These five loci yielded clear PCR products and were hence chosen over the remaining anonymous loci. One of our nine loci (02401) was chosen because it amplified best among a panel of 11 loci selected from a list of 242 anonymous markers proposed by Backström et al.  for population-genetic studies in birds. Another three loci, mitochondrial NADH dehydrogenase subunit 2 (ND2), β-fibrinogen intron 5 (Fib5) and tyrosinase related protein 1 (Tyrp1), are more conventional population-genetic and phylogenetic markers, two of which have amplified Elaenia DNA well in the past (ND2, Fib5; [20, 30]). Recent genomic data reveal that Tyrp1 is sex linked (present in the Z chromosome) in birds so sequence information of this locus was analyzed accordingly . General PCR conditions largely followed Rheindt et al. , with primers and annealing temperatures for the nine loci provided in Additional file 1: Table S4. We manually edited all sequence chromatograms following Rheindt et al. . Heterozygous sites were scored as degenerate sites based on established conventions of the International Union of Pure and Applied Chemistry (IUPAC), and were retained as one sequence per individual. All sequences generated in this study are available on GenBank (KY448474 - KY448797). We calculated the number of variable sites and parsimony-informative sites for our dataset in MEGA6 .
Tests for intragenic recombination and neutrality
We performed genetic analysis for recombinant detection (GARD)  in order to assess the presence of intragenic recombination for the nuclear loci within our dataset. We conducted separate GARD analyses for each locus in the Datamonkey webserver . We also performed the Hudson-Kreitman-Aguadé (HKA) test  as implemented in the program HKA https://bio.cst.temple.edu/~hey/software to assess deviations from neutrality for all loci in our dataset. We used the program SITES https://bio.cst.temple.edu/~hey/software  to generate input for the HKA test. The HKA test was carried out for the dataset consisting of the Elaenia taxa chilensis, albiceps and pallatangae only as this was our focal group for testing gene flow and isolation.
We calculated net between-taxa mean distances in MEGA using raw p-distances. We computed distances separately for the mitochondrial gene ND2 and for the all the autosomal nuclear loci, allowing for pairwise deletion of missing data. For the nuclear dataset, we computed genetic distances for each nuclear locus and averaged across all loci.
We reconstructed the phylogeny of the ‘montane Elaenia clade’ using both sequence concatenation as well as species tree methods . We concatenated sequences of all nuclear loci using Sequence Matrix v1.7.8  and employed RAxML GUI v1.5  to reconstruct a maximum likelihood (ML) tree under the GTR + gamma model, performing 100 runs using a bootstrap algorithm with 1000 replicates. The resulting tree was visualized in FigTree v1.4.2 . We performed species tree reconstructions using the maximum pseudo-likelihood coalescent method MP-EST  implemented in the STRAW  webserver, as well as a Bayesian method implemented in *BEAST v1.8.2 . The anonymous locus q25 could not be sequenced for E. fallax and hence was removed from species tree analyses.
For the MP-EST analysis we first obtained individual gene trees with 1000 bootstraps for the eight loci using RAxML following the same parameters as for the concatenated analysis. We further rooted each gene tree with the E. chiriquensis outgroup using the ReRoot module in the STRAW webserver. The species tree was then estimated in MP-EST. In the Bayesian species tree reconstruction we considered each locus as a separate partition. For each locus, we determined the substitution model through jModelTest v 2.1.7 [48, 49] (Additional file 1: Table S5). We used a separate relaxed molecular clock model for each locus and estimated relative clock rates, applying a Yule speciation process with random starting gene trees under two independent runs of 100 million generations each and sampling every 100th generation. We checked parameter convergence of both runs in Tracer v 1.6  and combined them in LogCombiner 1.8.2 with a 50% burnin while resampling once per 1000 generations. Subsequently we constructed a maximum clade credibility tree using median heights in TreeAnnotator 1.8.2 while allowing for a further 10% burnin using a posterior probability limit of 0.95. The final tree was visualized in FigTree. In order to explore reticulations and possible species tree topologies, we also used the combined trees obtained from LogCombiner 1.8.2 to generate a representation of topological uncertainty of species tree in DensiTree  with a 10% burnin.
Genetic clustering and admixture
We investigated gene flow among the three taxa of the ‘albiceps complex’ within the ‘montane Elaenia clade’ (i.e., chilensis, albiceps, pallatangae) with Bayesian clustering and admixture analyses in BAPS 5.4 [54–56]. We performed separate analyses on the alignment of the autosomal nuclear loci and the mitochondrial gene ND2. We first performed mixture analyses with ten iterations for a specified genetic cluster (K) and obtained estimates of ancestry coefficients of each individual. Based on our preliminary observations from phylogenetic analyses we performed mixture analyses for both K = 2 and K = 3. We further used the results from mixture analysis to perform population assignment and admixture analysis at both K = 2 and K = 3 using default settings. BAPS requires phased data for nuclear loci, hence we phased our data using PHASE v2.1 [57, 58] with default settings as implemented in DnaSP 5 .
Modeling gene flow with coalescent simulations
To further investigate gene flow among the three taxa chilensis, pallatangae and albiceps, we performed coalescent simulations of various evolutionary models in fastsimcoal2 v 188.8.131.52.21  and selected the best model using an Approximate Bayesian Computation (ABC) framework in the R package ABC . Fastsimcoal simulates data following the Wright-Fisher model of evolution assuming neutrality of the genetic markers, no recombination within loci and free recombination between loci and random mating. ABC based approaches can efficiently differentiate between genetic patterns arising from gene flow versus shared ancestral polymorphism [62–66]. Coalescent simulations were only carried out on loci that showed no significant evidence of intragenic recombination and deviation from neutrality. We simulated a no-gene flow model and various gene flow models in a three-population framework and used alternative topological arrangements resulting from species tree analyses as well as leveraging hypothesized earth-historic events to construct evolutionary models reflecting these scenarios.
We also removed the sex-linked locus Tyrp1 for this analysis and from the remaining set of eight loci (seven nuclear autosomal and one mitochondrial) we simulated sequences of the same length as obtained from our sequencing and alignments. We allowed for complete recombination between loci but no intra-locus recombination. We performed 1,000,000 simulations for each model and calculated summary statistics from both the observed data as well as simulated datasets using arlsumstats v. 3.5.2 . We performed separate simulations for the alignment of nuclear loci and the mitochondrial gene ND2. We set priors for population size (10 to 100,000) based on information obtained from Rheindt et al. . Following the species tree topology and net interspecies mitochondrial distances observed in this study, and accounting for the uncertainty of divergence estimates between these lineages, we applied broad priors for divergence time parameters (divergence between chilensis and pallatangae: 10,000 to 2,000,000; divergence between albiceps and ancestor of chilensis and pallatangae: 1,000,000 to 7,000,000). We chose a broad prior for migration based on our understanding of species biology (probability of migration: 1e-8 to 1e-4 per generation), with all priors following a uniform distribution. We performed all simulations assuming a constant population size. In our first set of models we assumed various continuous, symmetric and low gene flow scenarios alongside a no-gene flow model. We plotted prior and posterior distributions and verified that the simulations had effectively sampled from the prior distributions for our parameters of interest. We always considered chilensis and pallatangae as sister clades with a much more recent coalescence than the timing of coalescence between albiceps and the ancestral population of chilensis-pallatangae. Whenever necessary, we further refined our models to test more complex scenarios of asymmetric gene flow as well as of isolation and migration patterns.
We chose summary statistics that can explain differentiation between populations and their shared variability (Additional file 1: Table S6). Before performing model selection, we quantified model misclassification using both hard (confusion matrix; see ) and soft classification (mean posterior probability) schemes as implemented in the ABC package . We estimated the posterior probability of each model from a subset of simulations closest to the observed data (tolerance level) using multinomial logistic regression in the ABC package .
Ancestral state reconstruction
We performed ancestral character state reconstruction of underparts coloration using BayesTraits v 2.0 . We used a discrete two state coding (yellow versus white) for the terminal species nodes on the tree. The traits for each taxon are: pallatangae (yellow), chilensis (white), albiceps (white), frantzii (white), olivina (yellow), fallax (white), cherriei (white), martinica (white), chiriquensis (white). We used the multistate module in BayesTraits and performed a Bayesian MCMC run with 2,000,000 iterations, including a 100,000 iteration burnin, sampling every 500 generations, and estimated the ancestral state of each node in the final Bayesian species tree. We used an exponential distribution with a hyperprior as suggested by the program manual. We checked for parameter convergence using the coda package  in R.
Bioacoustic variation within the ‘albiceps complex’
We obtained 5573 bp of sequence data from across seven autosomal nuclear loci (3711 bp), one Z-linked locus (Tyrp1, 744 bp) and one mitochondrial locus (ND2, 1118 bp). Although we did not manage to sequence all 46 samples for all loci, we generated sequences for 67.4% of all sample–locus combinations (Additional file 1: Table S7). The number of variable sites varied between 13 and 246 per locus and the observed number of parsimony-informative sites per locus varied between 7 and 176 (Additional file 1: Table S8).
Recombination and neutrality
None of the loci showed any significant evidence of recombination. We performed pairwise HKA test comparisons and did not observe significant deviation from neutrality in any locus (all pairwise P values >0.01; Additional file 1: Table S9). Thus we proceeded to carry out coalescent simulations on all eight loci (total 4829 bp).
Raw mean pairwise p-distances for all nuclear (nDNA; below diagonal) and mitochondrial (mtDNA; above diagonal) DNA sequences
Bayesian clustering and gene flow
We truncated missing data at both terminal sequence ends and used the resulting alignment of 2802 bp of nuclear data for phasing to perform BAPS analysis. However, we did not use the phased alignment for the aforementioned phylogenetic reconstructions because the 909 bp loss in alignment length led to lower resolution. Clustering analyses of both nuclear and mitochondrial markers across the three members of the ‘albiceps complex’ revealed the presence of two major genetic clusters, both at K = 2 and K = 3, with chilensis and pallatangae belonging to a single cluster (Fig. 4b-e). Increasing cluster numbers from K = 2 to K = 3 did not add further obvious taxonomic or geographic resolution among members of the ‘albiceps complex’ (Fig. 4b-e).
Modeling gene flow
Posterior probabilities for all models based on multinomial logistic regression; Models tested using both nuclear loci and mitochondrial ND2
Model set 1
No gene flow
Gene flow between chilensis and pallatangae only
Gene flow between chilensis and albiceps only
Gene flow between both chilensis and pallatangae, and chilensis and albiceps
Gene flow among all three lineages
Posterior probabilities for all models based on multinomial logistic regression; models specifically tested for nuclear loci only
Model set 2
Equal gene flow between both chilensis and pallatangae, and chilensis and albiceps
Higher rate of gene flow between chilensis and pallatangae compared to chilensis and albiceps (asymmetric gene flow model)
Asymmetric gene flow model with gene flow until the last 20,000 years and isolation thereafter (isolation-migration model).
Ne albiceps (number of individuals)
Ne chilensis (number of individuals)
Ne pallatangae (number of individuals)
Time of divergence between chilensis and pallatangae (number of generations)
Time of divergence between chilensis pallatangae clade and albiceps (number of generations)
Migration rate between chilensis and pallatangae (probability of migration per generation)
Migration rate between chilensis and albiceps (probability of migration per generation)
1.24E + 04
2.41E + 04
4.24E + 02
5.56E + 04
1.11E + 06
3.12E + 04
4.92E + 04
3.11E + 04
1.62E + 05
2.62E + 06
7.66E + 04
8.50E + 04
7.38E + 04
6.69E + 05
5.78E + 06
7.35E + 04
8.20E + 04
7.14E + 04
7.11E + 05
5.52E + 06
9.47E + 04
9.50E + 04
7.97E + 04
4.68E + 05
6.67E + 06
9.91E + 04
9.93E + 04
9.82E + 04
1.51E + 06
6.95E + 06
1.00E + 05
1.00E + 05
1.00E + 05
1.98E + 06
7.00E + 06
Mitochondrial gene ND2
Ne albiceps (number of individuals)
Ne chilensis (number of individuals)
Ne pallatangae (number of individuals)
Time of divergence between chilensis and pallatangae (number of generations)
Time of divergence between chilensis pallatangae clade and albiceps (number of generations)
Migration rate between chilensis and pallatangae (probability of migration per generation)
6.83E + 01
3.97E + 01
1.12E + 01
1.11E + 04
1.00E + 06
3.57E + 03
4.80E + 03
3.67E + 03
1.24E + 05
1.10E + 06
1.90E + 04
1.91E + 04
1.86E + 04
9.32E + 05
2.55E + 06
1.75E + 04
1.79E + 04
1.72E + 04
9.65E + 05
2.99E + 06
2.31E + 04
2.34E + 04
2.33E + 04
7.31E + 05
1.77E + 06
2.47E + 04
2.48E + 04
2.47E + 04
1.93E + 06
6.56E + 06
2.50E + 04
2.50E + 04
2.50E + 04
2.00E + 06
6.99E + 06
Reconstruction of ancestral character states
The character state reconstruction analysis revealed that white underparts are a more probable ancestral character of the ‘albiceps complex’ (Fig. 3c). In general, for all internal nodes within the ‘montane Elaenia clade’, white underparts are more probable ancestral state than yellow underparts (Fig. 3c).
Recurrent Quaternary episodes of global cooling, the so-called “ice ages”, are known to have had an immense impact on differentiation processes at higher latitudes [4–8], yet very little is known about their effects on the biota of the tropics, including such centers of global biodiversity as the Andes. Using phylogenetic and population-genetic methods in conjunction with ~5.5 kb of sequence data across dozens of individuals of an Andean flycatcher lineage, coupled with morphological and bioacoustic data, our results suggest bouts of gene flow connecting some flycatcher lineages during peaks of global cooling when montane forest habitat would have extended to lower elevations and linked different mountain ridges [27, 28], while the current interglacial period of increased temperatures seems to be characterized by impeded gene flow between those lineages (Figs. 4 and 5; Tables 2 and 3).
Incongruence between species tree reconstructions and concatenation-based methods
Our concatenation-based tree of the montane Elaenia clade using eight nuclear loci and one mitochondrial gene (Fig. 3a) was in broad agreement with an earlier phylogenetic reconstruction based on only two loci that had also been concatenated , revealing a tight-knit grouping consisting of the forms pallatangae and chilensis, which, in turn, were deeply differentiated from the remaining six species-level lineages of the montane Elaenia clade. In contrast to the concatenation-based tree (Fig. 3a), one of our two species tree reconstructions (the one based on *BEAST; Fig. 3c) provides strong support for a placement of E. albiceps as sister to the pallatangae-chilensis duo, while the other species tree method (based on MP-EST) recovers the same arrangement but with negligible support (Fig. 3b). In these species tree arrangements, the three forms pallatangae, chilensis and albiceps together form what has been called the ‘albiceps complex’, a traditional taxonomic treatment that has long associated these three forms as close relatives (e.g. ). It is difficult to identify the causes for the strong disagreement between concatenation and species tree methods in the placement of E. albiceps. On the one hand, inherent methodological differences in these two approaches are known to lead to a frequent incidence of conflict between their resulting topologies [41, 72–75]. On the other hand, both approaches have different susceptibilities to methodological biases in terms of number of loci, average locus length, phasing and other factors [41, 76, 77].
In complicated radiations such as the ‘albiceps complex’ with a potentially rich history of inter-species gene flow, the future use of genome-wide DNA data is called for to achieve a final resolution of evolutionary trajectories. At the same time, our analyses allow for the conclusion that E. albiceps is a strongly diverged species distinct from E. pallatangae, and that the latter should include E. p. chilensis as a subspecies.
Phenotypic differentiation across the Andes in the face of gene flow
Species tree methods agreed with population-genetic approaches (e.g. mitochondrial network, BAPS) in aligning chilensis with pallatangae on a shallow node at near-zero mitochondrial differentiation (Figs. 3b, c and 4; Table 1). Contrary to Rheindt et al.’s  study of a single nuclear locus, our analysis of seven nuclear autosomal loci points to shallow levels of nuclear differentiation between chilensis and pallatangae (e.g. Table 1). In fact, our ABC analyses suggest considerable levels of recent gene flow between chilensis and pallatangae, at least up to a few thousand years ago (Tables 4 and 5). The two forms have never been considered conspecific because of their diverging underparts coloration (yellow versus white) providing them with superficial plumage differences that appear comparatively pronounced within the genus Elaenia. However, plumage differences are known to be less important in the reproductive isolation of tyrant-flycatcher species (Tyrannidae), especially as compared to bioacoustic differences . Our bioacoustic analysis showed that the calls of chilensis and pallatangae resemble each other more closely than either of them resembles the call note of E. albiceps, with clinal vocal trends in the pallatangae-chilensis cluster from north to south along the Andes (Fig. 2). Together, the genetic and phenotypic evidence leads us to propose that pallatangae and chilensis are two well-differentiated subspecies of a single biological species, E. pallatangae, connected by frequent gene flow and contiguity in vocal traits.
Gene flow during periods of glacial maxima
Our approximate Bayesian computations suggest a complicated evolutionary history of the ‘albiceps complex’ with variable bouts of gene flow among constituent members (model G, Fig. 5; Tables 2, 3, 4 and 5). While the detection of extensive gene flow between pallatangae and chilensis, along with their similar vocalizations, leads us to call for their inclusion within one biological species, there is also a genetic signal for reciprocal monophyly and much less regular gene flow between the taxa albiceps and chilensis (model G, Fig. 5; Tables 2, 3, 4 and 5). These two taxa are possibly parapatric in their breeding distributions, with albiceps in the northern and central Andes and chilensis in the southern Andes, but gene flow is conceivable in central Bolivia where they probably meet .
Paleoclimatological studies specifically from the La Paz and Cochabamba in Bolivia reveal that ice sheets had expanded greatly during the last ice age, pushing the tree line down the slope by at least 1000 m from its present elevation [27, 80–82]. Hence this region and Andean areas beyond would have provided habitat connectivity for various Elaenia species during glacial periods, thereby facilitating gene flow.
We conducted computations to test whether gene flow between chilensis and the other two forms of Elaenia flycatcher happens roughly uniformly across time (model F, Fig. 5; Tables 2 and 3), or whether there is a pronounced signal for gene flow to occur during periods of global cooling (=‘ice ages’) when Andean habitats shift down the slope in response to a drop in global temperatures (model G, Fig. 5; Tables 2 and 3; [27, 80–82]). The Quaternary has gone through multiple ice ages [1, 2, 27] that may have allowed for increased gene flow among Andean biota, but testing comparative coalescent models for all these ice ages is not possible for a dataset of our size as it may lack power to account for such complex scenarios. Therefore, we concentrated on the most recent period of global cooling (climaxing roughly at 20,000 years ago) that is most immediately previous to the present interglacial period of warmer global temperatures. Specifically, we tested whether the present interglacial would have been characterized by a lack of gene flow compared to bouts of gene flow during the previous period of global cooling. Our tests confirm that bouts of gene flow would have been more intense prior to the present interglacial period of warmer global temperatures (Table 3), attesting to the importance of Pleistocene montane habitat shifts in the Andes for gene flow among high-elevation biota.
While glacial periods are known to have been engines of population fragmentation at the colder latitudes of the northern hemisphere , this study is one of the first to demonstrate that glaciations may have had an opposite effect on super-diverse Neotropical mountain ranges, providing montane habitat connections between populations on adjacent ridges that would otherwise be isolated from one another for millions of years. On the surface, this mechanism seems to counteract speciation by providing recurrent gene flow between incipient lineages, but it may actually be one of the underlying reasons for the unparalleled species diversity of the Andes: Montane forest bridges during ice ages may help maintain the exceptional species diversity of the Andes and other tropical mountain ranges, while ice ages would expose temperate mountain chains to the spread of less hospitable cool-adapted habitats, such as treeline scrub, tundra or ice sheets, that would lead to widespread extinction events of forest biota.
Phenotypic character displacement
Our analyses uncover a complicated pattern of gene flow between the taxon chilensis and two other closely related forms of Elaenia flycatcher. While the white-bellied chilensis is connected with the yellow-bellied pallatangae by extensive gene flow and vocal similarities (Figs. 2 and 4; Tables 2, 3, 4 and 5; Additional file 1: Tables S10 and S12) – to the extent that we place them in the same biological species – there seems to be sporadic gene flow also between white-bellied chilensis and white-bellied albiceps. Vocal differences between albiceps and the other two forms are more pronounced (Fig. 2; Table 2, Additional file 1: Tables S10A and S12A). When one additionally considers the near-complete sympatry between albiceps and pallatangae (Fig. 1), it becomes clear that albiceps is an independent species, and rare instances of gene flow between it and any other Elaenia taxon must be treated as inter-specific introgression .
We studied the effects of Pleistocene climatic fluctuations on speciation dynamics in the ‘albiceps complex’ of Andean Elaenia tyrant-flycatchers. Using multi-locus sequence data as well as morphological and bioacoustics data we asked if ice ages have affected diversification within this group through habitat shifts in the montane biome. Our analyses reveal that periods of global cooling have bolstered gene flow among parapatric lineages. We also found strong evidence of reproductive isolation between sympatric lineages, and a pattern of character displacement in vocal and plumage traits from areas of sympatry to those of allopatry. Our findings provide a compelling mechanistic explanation for diversification and character displacement in natural populations within the Andean biodiversity hotspot.
We thank Miguel Alcaide, Maude Baldwin, Ingrid Soltero, Niclas Backström and Matt Fujita for help and advice in the laboratory. We would like to acknowledge the following people at the following institutions for the provision of tissue samples: Donna Dittmann, Van Remsen and Robb Brumfield (Louisiana State University Museum of Natural Science, Baton Rouge, Louisiana), Mark Robbins and A. Townsend Peterson (University of Kansas Natural History Museum, Lawrence, Kansas), David Willard and Shannon Hackett (Field Museum of Natural History, Chicago, Illinois), Paul Sweet (American Museum of Natural History, New York), Christopher Huddleston (Smithsonian Institution National Museum of Natural History, Washington, D.C.), and Jon Fjeldså (Zoological Museum, University of Copenhagen).
FER was funded by a Singapore Ministry of Education Tier I research grant (WBS R-154-000-658-112). BC acknowledges support by SEABIG grants (grant numbers: WBS R-154-000-648-646 and WBS R-154-000-648-733). We acknowledge a postdoctoral fellowship by the DAAD for having funded some of the work performed.
Availability of data and materials
All sequence data are available in GenBank (accession numbers: KY448474 - KY448797).
FER and SVE conceived the study and research questions, FER obtained the samples and generated all the sequence data, BC and KG performed all morphological and genetic analyses, CYG collected acoustic data and performed bioacoustic analyses. BC and FER drafted the manuscript with input from SVE, KG and CYG. All authors have approved of the final version of the manuscript.
Ethics approval and consent to participate
All samples used in this study were obtained from various Natural History museum collections (see Acknowledgments) under proper tissue loans and import permits.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Bintanja R, Van Oldenborgh G, Drijfhout S, Wouters B, Katsman C. Important role for ocean warming and increased ice-shelf melt in Antarctic sea-ice expansion. Nat Geosci. 2013;6(5):376–9.View ArticleGoogle Scholar
- Johnson NK, Cicero C. New mitochondrial DNA data affirm the importance of Pleistocene speciation in North American birds. Evolution. 2004;58(5):1122–30.View ArticlePubMedGoogle Scholar
- Weigelt P, Steinbauer MJ, Cabral JS, Kreft H. Late quaternary climate change shapes island biodiversity. Nature. 2016;532:99–102.View ArticlePubMedGoogle Scholar
- Alsos IG, Engelskjøn T, Gielly L, Taberlet P, Brochmann C. Impact of ice ages on circumpolar molecular diversity: insights from an ecological key species. Mol Ecol. 2005;14(9):2739–53.View ArticlePubMedGoogle Scholar
- Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405(6789):907–13.View ArticlePubMedGoogle Scholar
- Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF. Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998;7(4):453–64.View ArticlePubMedGoogle Scholar
- Kotlík P, Deffontaine V, Mascheretti S, Zima J, Michaux JR, Searle JB. A northern glacial refugium for bank voles (Clethrionomys glareolus). Proc Natl Acad Sci USA. 2006;103(40):14860–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Seddon J, Santucci F, Reeve N, Hewitt G. DNA footprints of European hedgehogs, Erinaceus europaeus and E. concolor: Pleistocene refugia, postglacial expansion and colonization routes. Mol Ecol. 2001;10(9):2187–98.View ArticlePubMedGoogle Scholar
- Bowie RC, Fjeldså J, Hackett SJ, Bates JM, Crowe TM. Coalescent models reveal the relative roles of ancestral polymorphism, vicariance, and dispersal in shaping phylogeographical structure of an African montane forest robin. Mol Phylogenet Evol. 2006;38(1):171–88.View ArticlePubMedGoogle Scholar
- Bowie RC, Fjeldså J, Hackett SJ, Crowe TM. Molecular evolution in space and through time: mtDNA phylogeography of the Olive Sunbird (Nectarinia olivacea/obscura) throughout continental Africa. Mol Phylogenet Evol. 2004;33(1):56–74.View ArticlePubMedGoogle Scholar
- Lessa EP, Cook JA, Patton JL. Genetic footprints of demographic expansion in North America, but not Amazonia, during the late quaternary. Proc Natl Acad Sci USA. 2003;100(18):10331–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Haffer J. Speciation in Amazonian forest birds. Science. 1969;165(3889):131–7.View ArticlePubMedGoogle Scholar
- Burns KJ, Naoki K. Molecular phylogenetics and biogeography of Neotropical tanagers in the genus Tangara. Mol Phylogenet Evol. 2004;32(3):838–54.View ArticlePubMedGoogle Scholar
- Brumfield RT, Edwards SV. Evolution into and out of the Andes: a Bayesian analysis of historical diversification in Thamnophilus antshrikes. Evolution. 2007;61(2):346–67.View ArticlePubMedGoogle Scholar
- Elmer KR, Bonett RM, Wake DB, Lougheed SC. Early Miocene origin and cryptic diversification of South American salamanders. BMC Evol Biol. 2013;13(1):1.View ArticleGoogle Scholar
- Rull V. Speciation timing and neotropical biodiversity: the tertiary–quaternary debate in the light of molecular phylogenetic evidence. Mol Ecol. 2008;17(11):2722–9.View ArticlePubMedGoogle Scholar
- Santos JC, Coloma LA, Summers K, Caldwell JP, Ree R, Cannatella DC. Amazonian amphibian diversity is primarily derived from late Miocene Andean lineages. PLoS Biol. 2009;7(3):e1000056.View ArticlePubMed CentralGoogle Scholar
- Voelker G, Outlaw RK, Bowie RC. Pliocene forest dynamics as a primary driver of African bird speciation. Glob Ecol Biogeogr. 2010;19(1):111–21.View ArticleGoogle Scholar
- Rheindt FE, Christidis L, Cabanne GS, Miyaki C, Norman JA. The timing of neotropical speciation dynamics: a reconstruction of Myiopagis flycatcher diversification using phylogenetic and paleogeographic data. Mol Phylogenet Evol. 2009;53(3):961–71.View ArticlePubMedGoogle Scholar
- Rheindt FE, Christidis L, Norman JA. Habitat shifts in the evolutionary history of a Neotropical flycatcher lineage from forest and open landscapes. BMC Evol Biol. 2008;8(1):1.View ArticleGoogle Scholar
- Roy MS. Recent diversification in African greenbuls (Pycnonotidae: Andropadus) supports a montane speciation model. Proc R Soc B. 1997;264(1386):1337–44.View ArticlePubMed CentralGoogle Scholar
- Weir JT. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution. 2006;60(4):842–55.View ArticlePubMedGoogle Scholar
- Moyle RG, Manthey JD, Hosner PA, Rahman M, Lakim M, Sheldon FH. A genome-wide assessment of stages of elevational parapatry in Bornean passerine birds reveals no introgression: implications for processes and patterns of speciation. PeerJ. 2017;5:e3335.View ArticlePubMedPubMed CentralGoogle Scholar
- Hoorn C, Wesselingh F, Ter Steege H, Bermudez M, Mora A, Sevink J, Sanmartín I, Sanchez-Meseguer A, Anderson C, Figueiredo J. Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity. Science. 2010;330(6006):927–31.View ArticlePubMedGoogle Scholar
- Myers N, Mittermeier RA, Mittermeier CG, Da Fonseca GA, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403(6772):853–8.View ArticlePubMedGoogle Scholar
- Jouzel J, Masson-Delmotte V, Cattani O, Dreyfus G, Falourd S, Hoffmann G, Minster B, Nouet J, Barnola J-M, Chappellaz J. Orbital and millennial Antarctic climate variability over the past 800,000 years. Science. 2007;317(5839):793–6.View ArticlePubMedGoogle Scholar
- Hooghiemstra H. Quaternary and Upper-Pliocene glaciations and forest development in the tropical Andes: evidence from a long high-resolution pollen record from the sedimentary basin of Bogotá, Colombia. Palaeogeogr Palaeoclimatol Palaeoecol. 1989;72:11–26.View ArticleGoogle Scholar
- Newton I. Speciation and biogeography of birds. London: Academic Press; 2003.Google Scholar
- Fjeldså J, Bowie RC, Rahbek C. The role of mountain ranges in the diversification of birds. Annu Rev Ecol Evol Syst. 2012;43:249–65.View ArticleGoogle Scholar
- Rheindt FE, Christidis L, Norman JA. Genetic introgression, incomplete lineage sorting and faulty taxonomy create multiple cases of polyphyly in a montane clade of tyrant-flycatchers (Elaenia, Tyrannidae). Zool Scr. 2009;38(2):143–53.View ArticleGoogle Scholar
- Ridgely RS, Tudor G. Field guide to the songbirds of South America: the passerines. Austin: University of Texas Press; 2009.Google Scholar
- R Core Team R. A language and environment for statistical computing. Vienna https://www.r-project.org/: R Foundation for Statistical Computing; 2016.Google Scholar
- Jennings WB, Edwards SV. Speciational history of Australian grass finches (Poephila) inferred from thirty gene trees. Evolution. 2005;59(9):2033–47.PubMedGoogle Scholar
- Backström N, Fagerberg S, Ellegren H. Genomics of natural bird populations: a gene-based set of reference markers evenly spread across the avian genome. Mol Ecol. 2008;17(4):964–80.View ArticlePubMedGoogle Scholar
- Nadeau NJ, Burke T, Mundy NI. Evolution of an avian pigmentation gene correlates with a measure of sexual selection. Proc R Soc B. 2007;274(1620):1807–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Pond SLK, Posada D, Gravenor MB, Woelk CH, Frost SD. GARD: a genetic algorithm for recombination detection. Bioinformatics. 2006;22(24):3096–8.View ArticleGoogle Scholar
- Delport W, Poon AF, Frost SD, Pond SLK. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics. 2010;26(19):2455–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Hudson RR, Kreitman M, Aguadé M. A test of neutral molecular evolution based on nucleotide data. Genetics. 1987;116(1):153–9.PubMedPubMed CentralGoogle Scholar
- Hey J and Wakeley J. A coalescent estimator of the population recombination rate. Genet. 1997;145:833–46.Google Scholar
- Edwards SV. Is a new and general theory of molecular systematics emerging? Evolution. 2009;63(1):1–19.View ArticlePubMedGoogle Scholar
- Vaidya G, Lohman DJ, Meier R. SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information. Cladistics. 2011;27(2):171–80.View ArticleGoogle Scholar
- Silvestro D, Michalak I. raxmlGUI: a graphical front-end for RAxML. Org Divers Evol. 2012;12(4):335–7.View ArticleGoogle Scholar
- Rambaut A. FigTree v1. 4.2: Tree figure drawing tool. 2014.Google Scholar
- Liu L, Yu L, Edwards SV. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evol Biol. 2010;10(1):302.View ArticlePubMedPubMed CentralGoogle Scholar
- Shaw TI, Ruan Z, Glenn TC, Liu L. STRAW: species TRee analysis web server. Nucleic Acids Res. 2013;41:W238–41.View ArticlePubMedPubMed CentralGoogle Scholar
- Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29(8):1969–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.View ArticlePubMedPubMed CentralGoogle Scholar
- Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52(5):696–704.View ArticlePubMedGoogle Scholar
- Rambaut A, Suchard M, Xie D, Drummond A. Tracer v1. 6. 2014.Google Scholar
- Bouckaert RR. DensiTree: making sense of sets of phylogenetic trees. Bioinformatics. 2010;26(10):1372–3.View ArticlePubMedGoogle Scholar
- Leigh JW, Bryant D. popart: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6(9):1110–6.View ArticleGoogle Scholar
- Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–9.View ArticlePubMedGoogle Scholar
- Corander J, Marttinen P. Bayesian identification of admixture events using multilocus molecular markers. Mol Ecol. 2006;15(10):2833–43.View ArticlePubMedGoogle Scholar
- Corander J, Marttinen P, Sirén J, Tang J. Enhanced Bayesian modelling in BAPS software for learning genetic structures of populations. BMC Bioinf. 2008;9(1):1.View ArticleGoogle Scholar
- Tang J, Hanage WP, Fraser C, Corander J. Identifying currents in the gene pool for bacterial populations using an integrative approach. PLoS Comput Biol. 2009;5(8):e1000455.View ArticlePubMedPubMed CentralGoogle Scholar
- Stephens M, Donnelly P. A comparison of bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet. 2003;73(5):1162–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68(4):978–89.View ArticlePubMedPubMed CentralGoogle Scholar
- Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.View ArticlePubMedGoogle Scholar
- Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M. Robust demographic inference from genomic and SNP data. PLoS Genet. 2013;9(10):e1003905.View ArticlePubMedPubMed CentralGoogle Scholar
- Csilléry K, François O, Blum MG. abc: an R package for approximate Bayesian computation (ABC). Methods Ecol Evol. 2012;3(3):475–9.View ArticleGoogle Scholar
- Raposo do Amaral F, Albers PK, Edwards SV, Miyaki CY. Multilocus tests of Pleistocene refugia and ancient divergence in a pair of Atlantic Forest antbirds (Myrmeciza). Mol Ecol. 2013;22(15):3996–4013.View ArticlePubMedGoogle Scholar
- Resende-Moreira LC, Ramos ACS, Scliar MO, Silva RM, Azevedo VC, Ciampi AY, Lemos-Filho JP, Lovato MB. Gene flow between vicariant tree species: insights into savanna-forest evolutionary relationships. Tree Genet Genomes. 2017;13(2):36.View ArticleGoogle Scholar
- Sousa V, Beaumont M, Fernandes P, Coelho M, Chikhi L. Population divergence with or without admixture: selecting models using an ABC approach. Heredity. 2012;108(5):521–30.View ArticlePubMedGoogle Scholar
- Turchetto C, Fagundes NJ, Segatto AL, Kuhlemeier C, Solis Neffa VG, Speranza PR, Bonatto SL, Freitas LB. Diversification in the South American Pampas: the genetic and morphological variation of the widespread Petunia axillaris complex (Solanaceae). Mol Ecol. 2014;23(2):374–89.View ArticlePubMedGoogle Scholar
- Zhou Y, Duvaux L, Ren G, Zhang L, Savolainen O, Liu J. Importance of incomplete lineage sorting and introgression in the origin of shared genetic variation between two closely related pines with overlapping distributions. Heredity. 2016;118:211–20.View ArticlePubMedPubMed CentralGoogle Scholar
- Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.View ArticlePubMedGoogle Scholar
- Pagel M, Meade A. BayesTraits v. 2.0. Reading: University of Reading; 2013.Google Scholar
- Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6(1):7–11.Google Scholar
- Rutter N, Coronato A, Helmens K, Rabassa J, Zárate M. Glaciations in North and South America from the Miocene to the Last Glacial maximum: comparisons, linkages and uncertainties. Dordrecht: Springer Science & Business Media; 2012.View ArticleGoogle Scholar
- Smith JA, Seltzer GO, Farber DL, Rodbell DT, Finkel RC. Early local last glacial maximum in the tropical Andes. Science. 2005;308(5722):678–81.View ArticlePubMedGoogle Scholar
- Degnan JH, Rosenberg NA. Discordance of species trees with their most likely gene trees. PLoS Genet. 2006;2(5):e68.View ArticlePubMedPubMed CentralGoogle Scholar
- Edwards SV, Liu L, Pearl DK. High-resolution species trees without concatenation. Proc Natl Acad Sci USA. 2007;104(14):5936–41.View ArticlePubMedPubMed CentralGoogle Scholar
- Kubatko LS, Degnan JH. Inconsistency of phylogenetic estimates from concatenated data under coalescence. Syst Biol. 2007;6(1):17–24.View ArticleGoogle Scholar
- Liu L, Yu L, Kubatko L, Pearl DK, Edwards SV. Coalescent methods for estimating phylogenetic trees. Mol Phylogenet Evol. 2009;53(1):320–8.View ArticlePubMedGoogle Scholar
- Camargo A, Avila LJ, Morando M, Sites JW Jr. Accuracy and precision of species trees: effects of locus, individual, and base pair sampling on inference of species trees in lizards of the Liolaemus darwinii group (Squamata, Liolaemidae). Syst Biol. 2011;61(2):272–88.View ArticlePubMedGoogle Scholar
- Corl A, Ellegren H. Sampling strategies for species trees: the effects on phylogenetic inference of the number of genes, number of individuals, and whether loci are mitochondrial, sex-linked, or autosomal. Mol Phylogenet Evol. 2013;67(2):358–66.View ArticlePubMedGoogle Scholar
- Rheindt FE, Norman JA, Christidis L. DNA evidence shows vocalizations to be a better indicator of taxonomic limits than plumage patterns in Zimmerius tyrant-flycatchers. Mol Phylogenet Evol. 2008;48(1):150–6.View ArticlePubMedGoogle Scholar
- Fitzpatrick J, Bates J, Bostwick K, Caballero I, Clock B, Farnsworth A, Hosner P, Joseph L, Langham G, Lebbin D. Family Tyrannidae (tyrant-flycatchers). Handb Birds World. 2004;9:170–462.Google Scholar
- Baker PA, Fritz SC. Nature and causes of quaternary climate variation of tropical South America. Quat Sci Rev. 2015;124:31–47.View ArticleGoogle Scholar
- Kull C, Imhof S, Grosjean M, Zech R, Veit H. Late Pleistocene glaciation in the Central Andes: temperature versus humidity control—a case study from the eastern Bolivian Andes (17 S) and regional synthesis. Glob Planet Chang. 2008;60(1):148–64.View ArticleGoogle Scholar
- Mark BG. Tracing tropical Andean glaciers over space and time: some lessons and transdisciplinary implications. Glob Planet Chang. 2008;60(1):101–14.View ArticleGoogle Scholar
- Rheindt FE, Edwards SV. Genetic introgression: an integral but neglected component of speciation in birds. Auk. 2011;128(4):620–32.View ArticleGoogle Scholar