Origin and cross-century dynamics of an avian hybrid zone
© The Author(s). 2017
Received: 12 April 2017
Accepted: 27 November 2017
Published: 15 December 2017
Characterizations of the dynamics of hybrid zones in space and time can give insights about traits and processes important in population divergence and speciation. We characterized a hybrid zone between tanagers in the genus Ramphocelus (Aves, Thraupidae) located in southwestern Colombia. We evaluated whether this hybrid zone originated as a result of secondary contact or of primary differentiation, and described its dynamics across time using spatial analyses of molecular, morphological, and coloration data in combination with paleodistribution modeling.
Models of potential historical distributions based on climatic data and genetic signatures of demographic expansion suggested that the hybrid zone likely originated following secondary contact between populations that expanded their ranges out of isolated areas in the Quaternary. Concordant patterns of variation in phenotypic characters across the hybrid zone and its narrow extent are suggestive of a tension zone, maintained by a balance between dispersal and selection against hybrids. Estimates of phenotypic cline parameters obtained using specimens collected over nearly a century revealed that, in recent decades, the zone appears to have moved to the east and to higher elevations, and may have become narrower. Genetic variation was not clearly structured along the hybrid zone, but comparisons between historical and contemporary specimens suggested that temporal changes in its genetic makeup may also have occurred.
Our data suggest that the hybrid zone likey resulted from secondary contact between populations. The observed changes in the hybrid zone may be a result of sexual selection, asymmetric gene flow, or environmental change.
Characterizations of hybrid zones allow one to make inferences about traits and processes relevant to understanding the origin and maintenance of differences between populations and species [1, 2]. A classic question about hybrid zones is how are they formed, with previous studies proposing two main hypotheses (reviewed by ). The hypothesis of secondary contact posits that hybrid zones result from expansion of populations that were previously isolated geographically and which interbreed in contact zones because complete reproductive isolation between them was not reached during the allopatric phase . An alternative hypothesis postulates that hybrid zones form in parapatry, by primary differentiation across ecological gradients . Secondary contact is likely if environments that presently allow the distributions of hybridizing populations to overlap were disjunct in the past, a scenario that predicts one should observe genetic signatures of demographic expansions. Alternatively, primary differentiation along a gradient would occur if the extent of suitable environments for the hybridizing populations has been stable over time; this predicts that populations have not expanded their ranges historically, and that the position of the hybrid zone (as indicated by the position of clines in molecular and morphological traits) is coincident with an environmental transition [1, 6].
Inferring hybrid-zone origins from current patterns of variation is challenging because, with time, genetic signatures of secondary contact or primary intergradation tend to erode [3, 7]. Alternatively, then, tests of hypotheses posed to account for the origin of hybrid zones may be conducted by examining the historical distribution of hybridizing taxa using paleodistributional modeling, an approach employing niche models that characterize the current distribution of species in climatic space to infer historical potential distributions given climatic conditions of the past [8, 9]. Such models of historical distributions represent hypotheses one can further test using molecular data to evaluate their population-genetic predictions, such as signatures of population growth for presumably expanding populations and of constant population size and isolation by distance in populations occurring within climatically stable areas [10, 11]. This approach has revealed that several hybrid zones likely originated following range expansions leading to secondary contact [12–15].
Another focus of studies on hybrid zones is the analysis of their temporal dynamics, which can allow understanding the role played by different evolutionary forces in such scenarios. When hybrid genotypes are less fit than parental genotypes, ‘tension zones’ are formed, which are maintained by a balance between the homogenizing effect of dispersal into the hybrid zone and the diversifying effect of selection against hybrids . If there is endogenous selection against hybrids, then there should be coincidence in location and concordance in width of clines describing the variation in different traits and loci across a hybrid zone, and such clines should remain stable over time [16, 17]. However, hybrid zones are often temporally dynamic (i.e. they may shift in location or change in width) and because clines for different traits may change in different ways, one can make inferences about the action of particular processes (e.g., natural selection, sexual selection, competition, asymmetric hybridization, dominance drive) based on dynamics observed for different characters . For example, discordant patterns of plumage and mitochondrial DNA variation across a hybrid zone between Setophaga warblers, coupled with behavioral experiments showing aggressive superiority of males of one species over the other, indicate that movement of this zone has likely been driven by competition-mediated asymmetric hybridization [19–21]. Temporal changes in the makeup of hybrid zones may also reflect natural or human-mediated environmental changes [22, 23].
The R. flammigerus system is particularly well suited to studying the role of different evolutionary forces at work in hybrid zones owing to the existence of variation in characters with different modes of inheritance, and, presumably, under different forms of selection. Variation in rump coloration across the hybrid zone likely reflects variation in the concentration of a single carotenoid pigment and is influenced by the environment because carotenoids are obtained from the diet [25, 32, 33]. Furthermore, based on the strong sexual dichromatism in this species, and the likelihood that its mating system involves some degree of polygamy , plumage coloration is probably influenced by sexual selection. In contrast, morphometric variation [25, 29] likely has a strong genetic basis [35–37] and could be subject to natural selection [38–40]. The value of considering traits or loci with different modes of inheritance and under different selective pressures to understand evolutionary forces at work in hybrid zones is illustrated by studies showing (1) that clines for traits involved in courtship are displaced with respect to clines for presumably neutral traits or loci, suggesting a role for sexual selection driving introgression [41, 42]; (2) that sex-linked molecular markers introgress over shorter distances than autosomal markers, suggesting a role for sex chromosomes in reproductive isolation [43, 44]; or (3) that there is more limited introgression in organellar DNA than in nuclear genes, suggesting selection acts more strongly on hybrids of the heterogametic sex (Haldane’s rule; ).
Here, we sought to evaluate whether the Ramphocelus hybrid zone in southwestern Colombia originated as a result of secondary contact or of primary differentiation, and to examine the zone’s dynamics over nearly a century to make inferences about the action of different evolutionary processes. To accomplish this, we (1) reconstructed the biogeographic and demographic history of the hybridizing populations based on ecological niche modeling and coalescent analyses of mtDNA sequence data, (2) characterized genetic, morphological and plumage variation across the hybrid zone, and (3) compared spatial patterns of variation in morphometrics, plumage coloration, and genetic structure between specimens collected at different times to assess possible changes in the position and width of the hybrid zone.
We characterized the Ramphocelus hybrid zone historically by examining specimens collected from 1894 to 1986 in the ornithological collections of the American Museum of Natural History, the Cornell University Museum of Vertebrates, Universidad del Valle, and the Instituto de Ciencias Naturales at Universidad Nacional de Colombia. To describe current patterns of variation, we collected 73 new specimens in 2007-2010: 65 of them are from localities ranging across the hybrid zone over a distance of 140 km by road connecting the cities of Cali and Buenaventura in department Valle del Cauca ; the remaining eight are from localities distant from the hybrid zone in the departments of Antioquia (3), Risaralda (3), and Cauca (2). Study skins and tissue samples are deposited in the Museo de Historia Natural de la Universidad de los Andes (ANDES, Additional file 1: Table S1). All specimen localities (historical and current) were plotted and their position along a transect line that best adjusted to points (estimated using a linear regression between latitude and longitude) was recorded. To construct character clines, we recorded the perpendicular position of each specimen on the regression line (Fig. 2; Additional file 1: Table S1) and calculated the distance from the northwest extreme of our study transect on the Pacific coast to the position of each specimen along the line.
To model potential distributions of the study taxa, we used 343 georeferenced localities obtained from museum specimens and reliable field observations from Costa Rica, Panama, Colombia, Ecuador, and Peru (; GBIF Data Portal, C. Sánchez, pers. comm., our observations). We associated localities with GIS layers for 19 climate variables at a c. 1 km resolution developed for the present (WorldClim; ). With these data, we used a maximum entropy approach (Maxent; ) based on current climate layers to generate models of the ecological niche and potential distribution of flammigerus and icteronotus at present. We conducted analyses in which localities of both forms and presumed hybrids were considered together to build potential distribution models, and also built separate models using data for flammigerus-like and icteronotus-like specimens; intermediate individuals were excluded from the latter analyses. Model performance in predicting present distributions (evaluated using receiver-operating-characteristic curves; ) was satisfactory (see below), validating the use of this approach to infer potential distributions in the past.
We projected models based on current climate data onto historical climate surfaces for 6000 years ago and for the Last Glacial Maximum (LGM; 21,000 years ago) to determine whether the distributions of our study taxa were likely disjunct in the past as predicted by the secondary contact hypothesis, or have likely been continuous as predicted by the primary intergradation hypothesis. This approach requires assuming ecological niche conservatism and that climate represents a long-term stable constraint on potential distributions. Because ecological niche models are based only on climatic data, they tend to overpredict potential distributions into areas where the study species do not occur owing to historical limitations to dispersal (e.g., the Amazon region in R. flammigerus). To reduce overprediction, we cropped maps of current and historical distributions to include only areas within the Andes Ecoregion . Although cropping potential distributions to this region probably did not remove all areas of model overprediction, it allowed for a semi-quantitative comparison of potential distributions across different time periods by calculating the extent of presence areas within the ecoregion. Maxent produces a continuous output ranging from zero to one describing the probability of the species potentially being present at different sites. We considered a threshold of 10% omission to categorize pixels as suitable or unsuitable for each of the time periods.
Genetic characterization and demographic history
We analyzed variation in DNA sequences of the cytochrome b mitochondrial gene for 58 of the flammigerus/icteronotus individuals collected in Colombia from 2007 to 2010. In addition, we obtained sequences for three individuals from Ecuador and two from Panama, and combined our data with two sequences of flammigerus available in GenBank: one from Ecuador (accession U15719.1; ) and one from Panamá (FJ799882.1; ; Additional file 1: Table S1).
DNA was extracted from tissue samples or toepad samples taken from specimens using a Qiagen DNeasy Tissue Kit or a phenol-chloroform method . PCRs used primers H16064 and L14996  in 24 μl amplification reactions using the following conditions: 42 ng of DNA, 0.416 mM dNTPs, 0.5 mM of each primer, 1.042 units of 10X buffer with 1.56 mM MgCl2, 0.0246 units/ml AmpliTaq DNA polymerase), and 16.5 of sterile ddH2O. Reactions began with an initial denaturation at 94 °C for 2 min, followed by 34 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 1 min, with a final extension phase at 72 °C for 7 min. PCR products were purified with Affymetrix Exosap-IT and sequenced in both directions. Sequences were edited, assembled and aligned using Geneious Pro 3.6.1 (http://www.geneious.com). The mean length of these sequences was 988.8 bp (range 888-1008 bp).
We also analyzed mtDNA from historical toepad samples of 87 specimens collected by C. G. Sibley in 1956 and housed at the Cornell University Museum of Vertebrates . DNA extraction and amplification of these samples was carried out in a historical DNA lab, following protocols to reduce the odds of contamination . For these specimens, we amplified and sequenced c. 210 bp of the cytochrome b gene (mean = 209.9, range = 203-210 bp).
We examined genealogical relationships among haplotypes observed in flammigerus and icteronotus at present using a maximum-likelihood (ML) phylogenetic analysis employing the GAMMA model. Nodal support was estimated using 1000 bootstrap replicates in RAxML, run from the RAxML BlackBox Web-Server . We used as outgroups sequences of Ramphocelus carbo (AF310048.1; ) and R. passerinii (EF529965.1; ).
To assess potential changes in the genetic makeup of the hybrid zone over time, we examined population structure separately for the 1956 specimens and for our samples collected in 2007-2010 (hereafter 2010 specimens) employing procedures implemented in the program ARLEQUIN v3.5 . Based on each data set, we conducted analyses of molecular variance (AMOVAs). We divided our sampling transect in three sectors of equal length: sector 1, 0 – 44 km; sector 2, 45-89 km; and sector 3, 90-134 km (Fig. 2). Within each sector, we grouped individuals collected within 1 km from each other in a single locality. We calculated F-statistics to estimate differentiation among sectors (FCT), among localities within sectors (FSC), and among localities among sectors (FST). Given that we lacked extensive sampling within point localities that would allow us to examine changes in genetic structure at a fine scale, this analysis allowed us to examine whether there has been any change in the way in which genetic variation is grossly distributed within and among different parts of the hybrid zone; if spatial genetic structure has become eroded (e.g., if introgressive hybridization has led to genetic homogenization across the transect), then one would expect an increase in genetic variation existing within sectors and a decrease in that existing among sectors over time (i.e., higher FCT in the past than at present). To make data comparable across time periods, sequences for 2010 specimens were trimmed to match the 210 bp available for the 1956 specimens. As an additional way to visualize potential changes in genetic structure over time, we constructed median-joining haplotype networks for the 1956 and 2010 samples using the program PopART (http://popart.otago.ac.nz/index.shtml).
To determine whether populations of flammigerus and icteronotus have experienced demographic expansions or if these taxa have exhibited historically stable population sizes, we examined trends in effective population size through time using the Extended Bayesian Skyline Plot (EBSP) method implemented in Beast v1.7.4  using sequence data for the 2010 specimens. Demographic expansions are expected if the hybrid zone originated following secondary contact and stable population sizes are expected if the zone originated by primary differentiation. Analyses used the HKY + Γ model, which was selected as the best fit to the data according to the Bayesian Information Criterion (BIC) in JModelTest v 2.1.3 . We ran the analysis for 25,000,000 iterations of which the first 10% were discarded as burn-in; genealogies and model parameters were sampled every 10,000 iterations. For time calibration we assumed a lognormal relaxed clock and a cytochrome-b substitution rate of 2.08% divergence per million years . We used the mean of the distribution of population size as a prior (parameter “demographic.populationMean”) calculated from a “Coalescent: constant time” tree prior, run with the same parameters as above. Because this analysis assumes no genetic structure within the sample, we only considered populations located between Cali and Buenaventura (i.e., from the hybrid-zone transect). We conducted an analysis in which we included data for all specimens obtained across the transect and also separate analyses for flammigerus-like and icteronotus-like specimens. As with niche modelling analyses, intermediate individuals were not included in analyses conducted separately by taxon. Skyline plots were built in R  with code written by Valderrama et al. .
To characterize the hybrid zone phenotypically, we measured six morphological characters on museum specimens (n = 139 males, 83 females) with dial calipers to the nearest 0.1 mm: wing length (chord of unflattened wing from bend of wing to longest primary), exposed culmen, bill depth (at the base), bill width (at the base), tail length (from point of insertion of central rectrices to tip of longest rectrix), and tarsus length (from the joint of tarsometatarsus and tibiotarsus to the lateral edge of last undivided scute). To describe morphological variation, we reduced variation in these characters using a principal components analysis (PCA).
We characterized plumage coloration based on reflectance spectra from 400 to 700 nm measured on the rump of adult museum specimens (n = 144 males, 70 females) using an Ocean Optics USB4F00243 Spectrometer with the SpectraSuite software (Ocean Optics). Three color measurements were estimated for each reflectance spectrum based on segment classification analysis : brightness, an index of how much light is reflected from the sample relative to a white standard; chroma, the saturation of color; and hue, which relates to the wavelength of maximum slope. These measurements were calculated using R code written by Parra .
Phenotypic clines and temporal dynamics
Here, D and C are the Y values at the plateau’s extremes (i.e. character values at each end of the hybrid zone in our case), B is a coefficient denoting the steepness of the curve, and E corresponds to the distance along the X axis where 50% of the value in Y is observed (also denoted ED50; ). We used the ED50 value as an estimate of cline center, and estimated cline width as the difference between the ED10 and ED90 values. We used this criterion to estimate cline width because ED10 and ED90 define the values of phenotypes in the X axis beyond which there are no intermediate individuals in morphology and color. For example, only individuals with scarlet rump would be observed in distances in X above ED90, whereas only yellow-rumped individuals would be observed in distances below ED10. We calculated the above parameters using the drm and ll.4 functions implemented in the drc package for R .
Due to differences in sampling effort over the hybrid zone across time periods, we evaluated the sensitivity of our estimates of cline center and width to sampling effects and calculated confidence limits for these parameters based on a bootstrapping procedure. For each time period and for both morphology and plumage reflectance data, we generated a bootstrap distribution of cline center and width estimated from 1000 data sets constructed by keeping sample size constant while resampling individuals with replacement. Because in some cases sample size was small, bootstrapping resulted in some samples in which variation was not clinal or in which estimates of cline parameters were unrealistic given the geographic extent of the hybrid zone; therefore, we only retained bootstrap samples in which the estimated cline center took values between 0 and 140 km. We compared estimates of cline center and width across the tree time periods using analyses of variance (ANOVA) and post-hoc Tukey tests treating estimates obtained in bootstrap samples as replicates. Likely due to a low number of individuals with morphological data for the western extreme of the hybrid zone in 2010, bootstrap estimates of cline width for this time period were highly variable and unrealistic; therefore, we do not report confidence limits for this parameter and did not include it in the ANOVA.
Validation of Hill-function methods
To validate the use of Hill functions to describe phenotypic and genetic clines and to infer cline center and width parameters, we compared our estimates based on Hill functions to estimates obtained by widely used cline-fitting algorithms (HZAR and Analyze) on a published data set that has been subject to cline analyses, namely the molecular data of seven loci of the Manacus hybrid zone in Panama . We used the published allele frequency data  to estimate cline centers and widths using our approach and compared such estimates to published parameters . Consistent estimates of parameters across methods would validate our Hill-function approach as an alternative method to estimate hybrid zone cline in scenarios where traditional algorithms cannot be used.
A potential distribution model developed under current climatic conditions in Maxent accurately predicted the present-day distributions of flammigerus + icteronotus, with an area under the ROC-curve score of 0.987. Because this suggests that the assumption that climate limits distributions in these taxa is reasonable, we projected models onto past climatic conditions to estimate the potential historical distributions of flammigerus + icteronotus as well as those of each taxon separately at 6000 and 21,000 years ago.
Population genetic structure in historical and recent specimens
Among localities among sectors
Among populations within sectors
Reduction of morphometric variation using PCA resulted in a first component (PC1) describing body size in both males and females. In both sexes, variables loading most heavily on this axis (which accounted for 24.3% of the variation in males and 34.7% in females) were tail length and wing chord. Thus, in the following we use PC1 as a general measurement of body size. We did not consider other principal component axes (e.g., PC2, on which bill dimensions loaded heavily) in additional analyses because they did not vary gradually across the hybrid zone.
Cline parameters estimated for morphological and coloration data in historical and recent specimens
Results of post-hoc ANOVA comparing of clines center (c) and widths (w) among three time periods for morphological (PC1) and plumage coloration (chroma) data using the bootstrap distributions of parameters infered from Hill functions
Patterns of variation in color in space and time were similar to those observed for morphology. Of the three measurements of plumage coloration, chroma showed the clearest clinal pattern of variation, ranging from the yellow icteronotus to the redder flammigerus (Fig. 6). The Hill-function estimates of cline centers did not differ significantly between the pre-1911 (72.9 km) and 1956-1986 (73.0 km) periods, but the cline-center estimate for 2010 (76.3 km) was significantly different, suggesting a displacement of around 3 km in eastward direction (Tables 3 and 4, Figs. 6h and 7, Additional file 2: Figure S1).
Estimates of cline width were substantially more uncertain, but also appear to differ between the present and past. Cline width has declined across time, being wider in the pre-1911 (width 41.7 km) than in the 1956-1986 (22.0 km) and 2010 (33.1 km) specimens (Table 3, Fig. 7). Likewise, the estimated cline width in chroma at present was c. 6.7 km, but was wider in the past: 31.5 km in the pre-1911 specimens and 32.4 km 1956-1986 specimens (Table 3, Fig. 7).
Because sample sizes for females were much lower than those of males and because variation among female specimens across the hybrid zone was not clearly clinal (e.g., Additional file 3: Figure S2) we did not estimate cline parameters for female morphology and plumage measurements. Likewise, because measurements of plumage hue and brightness for males and females did not show clear clinal trends (data not shown), we did not attempt to estimate cline parameters for these traits.
Based on patterns of genetic variation, fossil pollen data, and ecological niche modeling, several studies in the north temperate zone indicate that the origin of hybrid zones can be explained as a result of population expansions from isolated refugia during the Quaternary [12, 13, 71]. Although a similar hypothesis was proposed to account for the origin of contact zones in tropical rainforest organisms [72, 73], research on the origin of hybrid zones in the Neotropical region has been relatively limited [74, 75]. Our niche models indicate that potential distributions of R. f. icteronotus and R. f. flammigerus were likely disjunct at the LGM (21,000 ya), but were potentially in contact by 6000 ya, likely following a marked expansion of the geographic distribution of R. f. flammigerus. This scenario is supported by the historical demography analysis based on mtDNA sequence data, which indicates that populations have likely experienced significant range expansions, a pattern that awaits confirmation with multilocus data (see below) that should allow for lower uncertainty in estimates of population genetic parameters. Despite the broad credibility intervals around estimates of population size through time (except in the case of flammigerus, which shows strong evidence for expansion), taken together with results of niche modelling, the pattern observed in skyline plots is consistent with a scenario in which forms likely diverged in isolation (icteronotus in the Pacific lowlands and mid elevations of the Cordillera Occidental and flammigerus in refugial areas of the Cauca Valley) and then met as distributions expanded, presumably tracking the influence of Pleistocene climate change on vegetation [76, 77]. A recent study also suggested that historical climatic changes likely promoted changes in the geographic distributions of lowland Neotropical birds which are presently separated by the Andes . Our data also suggest that the divergence between the hybridizing Ramphocelus populations likely occurred in the Pleistocene, as indicated by low levels of mtDNA divergence suggesting recent differentiation. However, because Quaternary climatic oscillations started well before the LGM , it is possible that distribution ranges became disjunct and reconnected repeatedly at various times throughout the Pleistocene.
In contrast to our proposed scenario suggesting the origin of the Ramphocelus hybrid zone may date to at least 6000 before present, Sibley  hypothesized that contact between flammigerus and icteronotus resulted from recent anthropogenic deforestation and expansion of crops creating scrub and second-growth habitats, which are favored by these tanagers over dense rain forest. Although our analyses suggest that climatic conditions were suitable for contact between these forms thousands of years prior to major human-caused alterations in the area, it is likely that anthropogenic activities have facilitated contact between them, possibly leading to an increased incidence of hybridization in recent times.
A recent study on a hybrid zone between Heliconius butterflies located in the same geographic region where we studied hybridization in Ramphocelus also provided evidence consistent with the hypothesis of origin via secondary contact . Because there are additional documented cases of hybridization in the same general area of southwestern Colombia (e.g., other Heliconius butterfiles , Oophaga poison frogs ), work on the history of the region is necessary to better understand the origin and maintenance of hybrid zones across taxa [12, 15].
Our analyses are consistent with Sibley’s  overall characterization of the Ramphocelus hybrid zone: there is clinal variation in coloration and body size, with males exhibiting clearer trends than females (see also ). With the caveat that uncertainty around parameter estimates is broad, two main additional insights are provided by our cline analyses. First, our data consistently indicate that for each period, clines for morphology and chroma are coincident (i.e., they have equal or very similar centers). Second, variation in morphological PC1 and chroma followed the same trend over time: for both traits, clines appear to have moved slightly to the east and to have become narrower from the past to the present.
The coincidence of cline centers for different characters and their apparent concordance in width in each of the three periods is consistent with a tension-zone model. This is further supported by the width estimated for character clines. Assuming that the hybrid zone originated at least 6000 ya according to climate-based models, that generation time in Ramphocelus is 1-2 years, and that dispersal distances per generation lie somewhere between 1 and 20 km, the expected width of clines under a neutral diffusion model (equations in [7, 83]) would be between c. 25 km and more than 500 km. This is much wider than what we estimated (Table 2), which suggests some form of selection is acting. Alternatively, narrow clines may be a result of a more recent origin of the hybrid zone than implied by climate data. However, even if one assumes that the hybrid zone is as young as proposed by Sibley , the estimated clines appear narrower than expected under neutral diffusion (c. 4-75 km). In sum, our results lead us to hypothesize that the Ramphocelus hybrid zone is likely a tension zone, maintained by a balance between dispersal and barriers to gene flow . This interpretation contrasts with Sibley’s  conclusions that “gene exchange appears to be unimpeded” and that there is “no evidence of selection against the hybrids”, which he reached based on clinal patterns of variation in coloration and body size, and on the high abundance of putative hybrids. Future studies should further evaluate our tension-zone hypothesis by examining survival and mating success of intermediate phenotypes resulting from hybridization relative to parental types. In addition, detailed characterizations of patterns of variation along other contact zones between flammigerus and icteronotus located north of our study region would be of interest to determine whether intrinsic barriers to gene exchange between these taxa may exist independently of geographic or ecological context .
That cline centers do not appear to be coincident across different periods for each of the characters suggests movement of the Ramphocelus hybrid zone, but this interpretation needs to be tempered because samples were not taken at exactly the same localities at each time period and because the centers we estimated in some cases had wide support limits. In addition, our results should be interpreted with care given that the inferred geographic displacement of cline centers has been slight, corresponding to only c. 3 km between consecutive time periods. Nonetheless, our field observations indicate that the current patterns of variation are, in fact, different from those described by Sibley . Specifically, we frequently observed yellow-rumped individuals at Salado and Queremal, where Sibley did not report any, and we also have occasional records of yellow-rumped individuals near Cali, the eastern extreme of the transect. Thus, our quantitative analyses and field observations suggest that the icteronotus phenotype (yellow rump and smaller body size) has indeed extended to the east. Several previous studies have also reported on moving tension zones [85–87] although others documented spatial stability [88, 89].
A possible explanation for the concurrent movement of clines for different characters in hybrid zones is competitive advantage of one phenotype mediated by aggression or by sexual selection . Thus, movement of the Ramphocelus hybrid zone may have been driven by competition or sexual selection favoring the icteronotus phenotype. We believe it is unlikely that aggressive superiority of icteronotus explains this pattern as shown in other studies of hybrid zones owing to its smaller body size, although we note that in a Manacus (Pipridae) hybrid zone in Panama, the smaller M. vitellinus is more aggressive than the larger M. candei . Thus, it would be of interest to study mating patterns in the field and to conduct mate-choice experiments to determine whether the icteronotus phenotype has increased reproductive success as a result of sexual selection via male-male dominance or female choice [91, 92].
If sexual selection is based on carotenoid plumage color, which presumably is strongly influenced by the environment , then it is somewhat puzzling that coloration seems to have moved across the hybrid zone in concert with morphometric variation, which presumably has high heritability and is not known to be involved in mate choice. It is possible, however, that the ability to obtain, accumulate and metabolize carotenoids has a heritable genetic basis  and that the fitness advantages it confers may be linked to genes involved in reproductive isolation . If this were the case in R. flammigerus, then it would represent a plausible explanation for the likely movement of coloration in concert with others traits.
An alternative explanation for zone movement was proposed by Sibley , who predicted that the icteronotus phenotype would introgress across the hybrid zone as a consequence of increased gene flow from coastal to interior populations resulting from larger populations sizes in the former. This could be studied in the future with multilocus estimates of effective population sizes and of the magnitude of gene flow in both directions . In addition, if deforestation has indeed resulted in increases in population size as hypothesized by Sibley , then movement of the hybrid zone may also partly reflect anthropogenic influences. Because hybrid zones tend to become entrapped in areas of low population densities acting as sinks for migration [1, 89], any changes in population size related to habitat modification may have partly facilitated the observed movement of the Ramphocelus hybrid zone.
We note that the inferred slight movement of the hybrid zone involves not only a west-east displacement, but also a shift in its center to higher elevations in the Andes. Because the elevational ranges of tropical birds may shift upslope in response to global warming , it is also possible that the movement of the hybrid zone is related to climatic change over the past few decades [18, 22, 97, 98], as recently documented for hybrid zones between woodpecker (Picidae) and chickadee (Paridae) species in North America [99, 100].
In contrast to multiple studies on hybridization in birds finding significant mtDNA divergence between populations located away from the center of hybrid zones and clinal variation in haplotype frequencies across them (e.g., [41, 45, 55, 74, 88, 101–105]), mtDNA variation was not geographically structured in our study system, a likely consequence of recent divergence of the hybridizing populations or of high levels of introgression. In addition, because the probability of attaining genetic differentiation between isolated populations (i.e., significant differences in allele frequencies or even reciprocal monophyly) is a function of time and effective population sizes , the lack of substantial mtDNA divergence along the Ramphocelus hybrid zone may have resulted from a short duration of the allopatric phase relative to effective population sizes. Because no clinal mtDNA variation was observed across the Ramphocelus hybrid zone, we were unable to estimate cline parameters for different time periods to examine zone movement as done in a few studies on hybrid zones examining genetic variation in specimens collected at different times [23, 87, 103, 107]. However, our analyses revealing more limited genetic structure across the zone as indicated by F-statistics in 2010 relative to 1956 provide some preliminary evidence that patterns of genetic variation may have not been stable over time. We hypothesize that these results may reflect an increase in introgression of mtDNA across our study transect since Sibley’s  time, but we acknowledge that because we only assayed a small fragment of mtDNA and because sampling was not spatially even between time periods, drawing any conclusion at this time would be premature. Nonetheless, our preliminary genetic data suggest that assessing temporal variation in spatial genetic structure using genome-wide markers (cf. ) would represent a fruitful avenue to better understand the ecological and evolutionary forces at work in this moving hybrid zone. Also, if shallow mtDNA divergence reflects overall patterns in genome-wide neutral divergence, then the Ramphocelus hybrid zone appears especially well-suited for analyses of variation in functionally important genes, which may provide important insights about the genetic basis of phenotypic variation. Recent work on other avian systems has revealed that phenotypic divergence in plumage traits controlled by relatively few loci may persist despite high overall genomic similarity [108–110].
Finally, our study serves to illustrate that Hill functions and bootstrap resampling are useful statistical tools to characterize patterns of variation across hybrid zones. In contrast to other studies , our work was not systematically designed to examine changes in the structure of the Ramphocelus hybrid zone over time, with repeated sampling of multiple specimens at fixed locations; rather, specimens were collected by numerous researchers, often opportunistically, for different purposes, and with no common sampling schemes at different times. Therefore, we did not feel confident in assigning individuals to discrete localities as required for cline-fitting algorithms which rely on calculating means and variances for samples of specimens from the same sites [66–68]. Our reanalysis of a previously published data set suggests that cline parameters estimated using Hill functions mirror closely those estimated by software commonly employed in studies of hybrid zones, suggesting such functions are useful alternatives when methods designed to study clines cannot be employed due to the nature of the available data. Although Hill functions have been used mainly to model dose-response curves, in principle they can be used to describe any relationship between variables that follows a sigmoid shape. In addition to not requiring data to be grouped in localities (i.e. each specimen conserves its position along sampling transects), the approach we followed does not assume normality in the data and can be run with low computational requeriments and relatively small sample sizes. On the other hand, we suspected that due to variation in the spatial distribution and number of specimens over time, inferences about temporal changes in the structure of the hybrid zone could be compromised by unequal sampling. However, by analyzing variation in bootstrap samples of specimens, we were able to account for sampling effects, finding that despite inconsistent sampling over time, there is sufficient signal in our data to document temporal changes in spatial patterns of phenotypic variation. Similar approaches could be employed to analyse temporal changes in spatial variation in genetic or phenotypic traits in other systems in which museum specimens have not been systematically collected over time and space but nonetheless contain rich historical information.
Models of potential historical distributions based on climatic data and genetic signatures of demographic expansion suggested that the Ramphocelus hybrid zone we studied originated following secondary contact between populations that expanded their ranges out of isolated areas in the Quaternary, likely as a consequence of climatic change. Concordant patterns of variation in phenotypic characters across the hybrid zone and its narrow extent are suggestive of a tension zone, maintained by a balance between dispersal and selection against hybrids. In addition, our estimates of phenotypic cline parameters obtained using specimens collected over nearly a century revealed that the zone has likely moved to the east and to higher elevations, and has likely become narrower. These observed changes in the hybrid zone may be a result of sexual selection, asymmetric gene flow, or environmental change. Our data represent a baseline for a variety of ecological, behavioral, and genetic studies that could shed further light on the forces involved in speciation and hybrid-zone dynamics in this system.
The Ministerio de Ambiente, Vivienda y Desarrollo Territorial of Colombia provided permits for research, collecting, and accessing genetics resources (RGE0014). F. Ayerbe, J. P. Lopéz, J. P.Goméz, N. Gutiérrez, S. Gonzalez, P. Tovar, F. Gonzalez, C. Wagner, J. Sandoval, the Luna Family, and the Castaño family helped in the field. For assistance with analyses, we are grateful to J. M. Cordovez, C. Ritz, C. Palacios, L. N. Céspedes, E. Derryberry, N. Gutiérrez, J. L. Parra, C. Pedraza, and E. Valderrama. We thank the Instituto de Génetica at Universidad de Los Andes (M. Linares) for access to facilities, and the curators and collection managers at institutions who allowed us to examine specimens or provided information or tissue samples for genetic analyses: American Museum of Natural History (J. Cracraft, P. Sweet), Cornell University Museum of Vertebrates (K. Botswick, I. Lovette), Instituto de Ciencias Naturales at Universidad Nacional de Colombia (F. G. Stiles), Louisiana State University Museum of Natural Science (R. Brumfield), Universidad del Cauca (M. D. P. Rivas), Universidad del Pacífico (A. Quintero), and Universidad del Valle (H. Alvarez-López).
The Facultad de Ciencias at Universidad de los Andes, an American Ornithologists’ Union Research Award, and a Collection Study Grant from the Department of Ornithology at the American Museum of Natural History provided funding for this study.
AMR and CDC conceived the study. AMR collected and measured specimens, conducted molecular genetic work, and performed initial analyses. MDC conducted molecular work on historical specimens. AMR and EAT analyzed data with supervision from CDC. AMR and CDC wrote the manuscript, with input from EAT and MDC. All authors read, commented on, and approved the final manuscript.
Ethics approval and consent to participate
Procedures employed have been approved by the ethics committee at Universidad de los Andes.
Consent for publication
All authors read, commented on, and approved the submission of the manuscript to BMC Evolutionary Biology.
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.
- Barton NH, Hewitt GM. Analysis of hybrid zones. Annu Rev Ecol Syst. 1985;16:113–48.View ArticleGoogle Scholar
- Coyne JA, Orr HA. Speciation. Sunderland: Sinauer Associates; 2004.Google Scholar
- Durrett R, Buttel L, Harrison R. Spatial models for hybrid zones. Heredity. 2000;84:9–19.PubMedView ArticleGoogle Scholar
- Mayr E. Systematics and the origin of species, from the viewpoint of a zoologist. Cambridge: Harvard University Press; 1942.Google Scholar
- Harrison RG. Hybrid zones: windows on evolutionary process. In: Futuyma DJ, Antonovics J, editors. Oxford surveys in evolutionary biology. Oxford: Oxford University Press; 1990. p. 69–128.Google Scholar
- Vines TH, Köhler SC, Thiel M, Ghira I, Sands TR, MacCallum CJ, et al. The maintenance of reproductive isolation in a mosaic hybrid zone between the fire-bellied toads Bombina bombina and B. variegata. Evolution. 2003;57:1876–88.PubMedView ArticleGoogle Scholar
- Endler JA. Geographic variation, speciation, and clines. Princeton: Princeton University Press; 1977.Google Scholar
- Peterson AT, Martínez-Meyer E, González-Salazar C. Reconstructing the Pleistocene geography of the Aphelocoma jays (Corvidae). Divers Distrib. 2004;10:237–46.View ArticleGoogle Scholar
- Carstens BC, Richards CL. Integrating coalescent and ecological niche modeling in comparative phylogeography. Evolution. 2007;61:1439–54.PubMedView ArticleGoogle Scholar
- Richards CL, Carstens BC, Knowles L. Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. J Biogeogr. 2007;34:1833–45.View ArticleGoogle Scholar
- Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C. Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot. Science. 2009;323:785–9.PubMedView ArticleGoogle Scholar
- Swenson NG. Gis-based niche models reveal unifying climatic mechanisms that maintain the location of avian hybrid zones in a North American suture zone. J Evol Biol. 2006;19:717–25.PubMedView ArticleGoogle Scholar
- Ruegg KC, Hijmans RJ, Moritz C. Climate change and the origin of migratory pathways in the Swainson’s thrush, Catharus ustulatus. J Biogeogr. 2006;33:1172–82.View ArticleGoogle Scholar
- Swenson NG. The past and future influence of geographic information systems on hybrid zone, phylogeographic and speciation research. J Evol Biol. 2008;21:421–34.PubMedView ArticleGoogle Scholar
- Moritz C, Hoskin CJ, MacKenzie JB, Phillips BL, Tonione M, Silva N, et al. Identification and dynamics of a cryptic suture zone in tropical rainforest. Proc R Soc Lond B Biol Sci. 2009;276:1235–44.View ArticleGoogle Scholar
- Moore BR, Price JT. Nature of selection in the northern flicker hybrid zone and its implications for speciation theory. In: Harrison RG, editor. Hybrid zones and the evolutionary process. Oxford: Oxford University Press; 1993. p. 196–225.Google Scholar
- Kruuk LEB, Baird SJE, Gale KS, Barton NH. A comparison of multilocus clines maintained by environmental adaptation or by selection against hybrids. Genetics. 1999;153:1959–71.PubMedPubMed CentralGoogle Scholar
- Buggs RJA. Empirical study of hybrid zone movement. Heredity. 2007;99:301–12.PubMedView ArticleGoogle Scholar
- Rohwer S, Bermingham E, Wood C. Plumage and mitochondrial DNA haplotype variation across a moving hybrid zone. Evolution. 2001;55:405–22.PubMedView ArticleGoogle Scholar
- Krosby M, Rohwer S. A 2000 km genetic wake yields evidence for northern glacial refugia and hybrid zone movement in a pair of songbirds. Proc R Soc Lond B Biol Sci. 2009;276:615–21.View ArticleGoogle Scholar
- Krosby M, Rohwer S. Ongoing movement of the Hermit Warbler X Townsend’s Warbler hybrid zone. PLoS One. 2010;5:e14164.PubMedPubMed CentralView ArticleGoogle Scholar
- Taylor SA, Larson EL, Harrison RG. Hybrid zones: windows on climate change. Trends Ecol Evol. 2015;30:398–406.PubMedPubMed CentralView ArticleGoogle Scholar
- Leaché AD, Grummer JA, Harris RB, Breckheimer IK. Evidence for concerted movement of nuclear and mitochondrial clines in a lizard hybrid zone. Mol Ecol. 2017;23:75.Google Scholar
- Griscom L. Notes on imaginary species of Ramphocelus. Auk. 1932;49:199–203.View ArticleGoogle Scholar
- Sibley CG. Hybridization in some colombian tanagers, avian genus Ramphocelus. Proc Am Philos Soc. 1958;102:448–53.Google Scholar
- Olson SL, Violani C. Some unusual hybrids of Ramphocelus, with remarks on evolution in the genus (Aves: Thraupinae). Boll Mus Regionale Sci Nat Torino. 1995;13:297–312.Google Scholar
- McCarthy EM. Handbook of avian hybrids of the world. Oxford: Oxford University Press, USA; 2006.Google Scholar
- Chapman FM. The distribution of bird-life in Colombia; a contribution to a biological survey of South America. Bull Am Mus Nat Hist. 1917;36:1–729.Google Scholar
- Bedoya MJ, Murillo OE. Evidencia morfológica de hibridación entre las subespecies de Ramphocelus flammigerus (Passeriformes: Thraupidae) en Colombia. Rev Biol Trop. 2012;60:75–85.PubMedView ArticleGoogle Scholar
- Remsen JV Jr, Areta JI, Cadena CD, Claramunt S, Jaramillo A, Pacheco JF, et al. A classification of the bird species of South America. American Ornithologists’ Union. 2017: Available from: http://www.museum.lsu.edu/~Remsen/SACCBaseline.htm.Google Scholar
- Hilty SL, Brown WL. A guide to the birds of Colombia. Princeton: Princeton University Press; 1986.Google Scholar
- Plath K. Color change in Ramphocelus flammigerus. Auk. 1945;62:304.Google Scholar
- Brush AH. Pigments in hybrid, variant and melanic tanagers (birds). Comp Biochem Physiol. 1970;36:785–93.View ArticleGoogle Scholar
- Krueger TR, Williams DA, Searcy WA. The genetic mating system of a tropical tanager. Condor. 2008;110:559–62.View ArticleGoogle Scholar
- Keller LF, Grant PR, Grant BR, Petren K. Heritability of morphological traits in Darwin’s finches: misidentified paternity and maternal effects. Heredity. 2001;87:325–36.PubMedView ArticleGoogle Scholar
- Bears H, Drever MC, Martin K. Comparative morphology of dark-eyed juncos Junco hyemalis breeding at two elevations: a common aviary experiment. J Avian Biol. 2008;39:152–62.View ArticleGoogle Scholar
- Ballentine B, Greenberg R. Common garden experiment reveals genetic control of phenotypic divergence between swamp sparrow subspecies that lack divergence in neutral genotypes. PLoS One. 2010;5:e10229.PubMedPubMed CentralView ArticleGoogle Scholar
- Schluter D, Smith JNM. Natural selection on beak and body size in the song sparrow. Evolution. 1986;40:221–31.PubMedView ArticleGoogle Scholar
- Grant PR, Grant BR. Unpredictable evolution in a 30-year study of Darwin’s finches. Science. 2002;296:707–11.PubMedView ArticleGoogle Scholar
- Benkman CW. Divergent selection drives the adaptive radiation of crossbills. Evolution. 2003;57:1176–81.PubMedView ArticleGoogle Scholar
- Brumfield RT, Jernigan RW, McDonald DB, Braun MJ. Evolutionary implications of divergent clines in an avian (Manacus: Aves) hybrid zone. Evolution. 2001;55:2070–87.PubMedView ArticleGoogle Scholar
- Baldassarre DT, White TA, Karubian J, Webster MS. Genomic and morphological analysis of a semipermeable avian hybrid zone suggests asymmetrical introgression of a sexual signal. Evolution. 2014;68:2644–57.PubMedView ArticleGoogle Scholar
- Carling MD, Brumfield RT. Speciation in Passerina buntings: introgression patterns of sex-linked loci identify a candidate gene region for reproductive isolation. Mol Ecol. 2009;18:834–47.PubMedView ArticleGoogle Scholar
- Taylor SA, Curry RL, White TA, Ferretti V, Lovette IJ. Spatiotemporally consistent genomic signatures of reproductive isolation in a moving hybrid zone. Evolution. 2014;68:3066–81.PubMedView ArticleGoogle Scholar
- Gowen FC, Maley JM, Cicero C, Peterson AT, Faircloth BC, Warr TC, et al. Speciation in Western Scrub-Jays, Haldane’s rule, and genetic clines in secondary contact. BMC Evol Biol. 2014;14:135.PubMedPubMed CentralView ArticleGoogle Scholar
- Darwin Database. BIOMAP alliance partners. 2007: Available from: http://www.biomap.net/BioMAP/login_en.php.Google Scholar
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.View ArticleGoogle Scholar
- Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.View ArticleGoogle Scholar
- Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ Conserv. 1997;24:38–49.View ArticleGoogle Scholar
- Dinerstein E. A conservation assessment of the terrestrial ecoregions of Latin America and the Caribbean. World Bank: Washington, D.C; 1995.View ArticleGoogle Scholar
- Hackett SJ. Molecular phylogenetics and biogeography of tanagers in the genus Ramphocelus (Aves). Mol Phylogenet Evol. 1996;5:368–82.PubMedView ArticleGoogle Scholar
- Burns KJ, Racicot RA. Molecular phylogenetics of a clade of lowland tanagers: implications for avian participation in the great American interchange. Auk. 2009;126:635–48.View ArticleGoogle Scholar
- Sambrook J, Russell DW. Molecular cloning: a laboratory manual. 3rd ed. Cold Spring Harbor: Cold Spring Harbor Laboratory; 2001.Google Scholar
- Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP. Primers for a PCR-based approach to mitochondrial genome sequencing in birds and other vertebrates. Mol Phylogenet Evol. 1999;12:105–14.PubMedView ArticleGoogle Scholar
- Carling MD, Serene LG, Lovette IJ. Using historical DNA to characterize hybridization between Baltimore Orioles (Icterus galbula) and Bullock’s Orioles (I. bullockii). Auk. 2011;128:61–8.View ArticleGoogle Scholar
- Stamatakis A, Hoover P, Rougemont J. A rapid bootstrap algorithm for the RAxML web servers. Syst Biol. 2008;57:758–71.PubMedView ArticleGoogle Scholar
- Sato A, Tichy H, O'hUigin C, Grant PR, Grant BR, Klein J. On the origin of Darwin’s finches. Mol Biol Evol. 2001;18:299–311.PubMedView ArticleGoogle Scholar
- Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.PubMedView ArticleGoogle Scholar
- Heled J, Drummond AJ. Bayesian inference of population size history from multiple loci. BMC Evol Biol. 2008;8:289.PubMedPubMed CentralView ArticleGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772.PubMedPubMed CentralView ArticleGoogle Scholar
- Weir JT, Price M. Andean uplift promotes lowland speciation through vicariance and dispersal in Dendrocincla woodcreepers. Mol Ecol. 2011;20:4550–63.PubMedView ArticleGoogle Scholar
- R Development Core Team. R: a language and environment for statistical computing. Vienna. URL http://www.r-project.org/: R Foundation for Statistical Computing; 2013.Google Scholar
- Valderrama E, Pérez-Emán JL, Brumfield RT, Cuervo AM, Cadena CD. The influence of the complex topography and dynamic history of the montane Neotropics on the evolutionary differentiation of a cloud forest bird (Premnoplex brunnescens, Furnariidae). J Biogeogr. 2014;41:1533–46.View ArticleGoogle Scholar
- Endler JA. On the measurement and classification of colour in studies of animal colour patterns. Biol J Linn Soc. 1990;41:315–52.View ArticleGoogle Scholar
- Parra JL. Color evolution in the hummingbird genus Coeligena. Evolution. 2009;64:324–35.PubMedView ArticleGoogle Scholar
- Derryberry EP, Derryberry GE, Maley JM, Brumfield RT. HZAR: hybrid zone analysis using an R software package. Mol Ecol Resour. 2014;14:652–63.PubMedView ArticleGoogle Scholar
- Barton NH, Baird S. Analyse: software for the analysis of geographic variation and hybrid zones. Edinburgh: University of Edinburgh; 1999.Google Scholar
- Porter AH, Wenger R, Geiger H, Scholl A, Shapiro AM. The Pontia daplidice-edusa hybrid zone in northwestern Italy. Evolution. 1997;51:1561–73.PubMedGoogle Scholar
- Ritz C, Baty F, Streibig JC, Gerhard D. Dose-response analysis using R. PLoS One. 2015;10:e0146021.PubMedPubMed CentralView ArticleGoogle Scholar
- Ritz C, Streibig JC. Bioassay analysis using R. J Stat Softw. 2005;12:1–22.View ArticleGoogle Scholar
- Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405:907–13.PubMedView ArticleGoogle Scholar
- Haffer J. Speciation in colombian forest birds west of the Andes. Am Mus Novit. 1967;2294:1–57.Google Scholar
- Endler JA. Problems in distinguishing historical from ecological factors in biogeography. Am Zool. 1982;22:441–52.View ArticleGoogle Scholar
- Brumfield RT. Mitochondrial variation in bolivian populations of the variable Antshrike (Thamnophilus caerulescens). Auk. 2005;122:414–32.View ArticleGoogle Scholar
- Dasmahapatra KK, Lamas G, Simpson F, Mallet J. The anatomy of a “suture zone” in Amazonian butterflies: a coalescent-based test for vicariant geographic divergence and speciation. Mol Ecol. 2010;19:4283–301.PubMedView ArticleGoogle Scholar
- Hooghiemstra H, Van der Hammen T. Quaternary ice-age dynamics in the Colombian Andes: developing an understanding of our legacy. Philos Trans R Soc Lond Biol Sci. 2004;359:173–81.View ArticleGoogle Scholar
- Ramírez-Barahona S, Eguiarte LE. The role of glacial cycles in promoting genetic diversity in the Neotropics: the case of cloud forests during the last glacial maximum. Ecol Evol. 2013;3:725–38.PubMedPubMed CentralView ArticleGoogle Scholar
- Cadena CD, Pedraza CA, Brumfield RT. Climate, habitat associations and the potential distributions of Neotropical birds: implications for diversification across the Andes. Rev Acad Colomb Cienc Ex Fis Nat. 2016;40:275.View ArticleGoogle Scholar
- Jansson R, Dynesius M. The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. Annu Rev Ecol Syst. 2002;33:741–77.View ArticleGoogle Scholar
- Arias CF, Rosales C, Salazar C, Castaño J, Bermingham E, Linares M, et al. Sharp genetic discontinuity across a unimodal Heliconius hybrid zone. Mol Ecol. 2012;21:5778–94.PubMedView ArticleGoogle Scholar
- Arias CF, Muñoz AG, Jiggins CD, Mavárez J, Bermingham E, Linares M. A hybrid zone provides evidence for incipient ecological speciation in Heliconius butterflies. Mol Ecol. 2008;17:4699–712.PubMedView ArticleGoogle Scholar
- Medina I, Wang IJ, Salazar C, Amézquita A. Hybridization promotes color polymorphism in the aposematic harlequin poison frog, Oophaga histrionica. Ecol Evol. 2013;3:4388–400.PubMedPubMed CentralView ArticleGoogle Scholar
- Barton NH, Gale KS. Genetic analysis of hybrid zones. In: Harrison RG, editor. Hybrid zones and the evolutionary process. New York: Oxford Univ. Press; 1993. p. 13–45.Google Scholar
- Harrison RG, Larson EL. Heterogeneous genome divergence, differential introgression, and the origin and structure of hybrid zones. Mol Ecol. 2016;25:2454–66.PubMedPubMed CentralView ArticleGoogle Scholar
- Gay L, Crochet P-A, Bell DA, Lenormand T. Comparing clines on molecular and phenotypic traits in hybrid zones: a window on tension zone models. Evolution. 2008;62:2789–806.PubMedView ArticleGoogle Scholar
- Dasmahapatra KK, Blum MJ, Aiello A, Hackwel S, Davies N, Bermingham EP, et al. Inferences from a rapidly moving hybrid zone. Evolution. 2002;56:741–53.PubMedView ArticleGoogle Scholar
- Smith KL, Hale JM, Gay L, Kearney M, Austin JJ, Parris KM, et al. Spatio-temporal changes in the structure of an Australian frog hybrid zone: a 40-year perspective. Evolution. 2013;67:3442–54.PubMedView ArticleGoogle Scholar
- Mettler RD, Spellman GM. A hybrid zone revisited: molecular and morphological analysis of the maintenance, movement, and evolution of a Great Plains avian (Cardinalidae: Pheucticus) hybrid zone. Mol Ecol. 2009;18:3256–67.PubMedPubMed CentralView ArticleGoogle Scholar
- Rosser N, Dasmahapatra KK, Mallet J. Stable Heliconius butterfly hybrid zones are correlated with a local rainfall peak at the edge of the Amazon basin. Evolution. 2014;68:3470–84.PubMedView ArticleGoogle Scholar
- McDonald DB, Clay RP, Brumfield RT, Braun MJ. Sexual selection on plumage and behavior in an avian hybrid zone: experimental tests of male-male interactions. Evolution. 2001;55:1443–51.PubMedView ArticleGoogle Scholar
- Bronson CL, Grubb TCJ, Sattler DG, Braun MJ. Mate preference: a possible causal mechanism for a moving hybrid zone. Anim Behav. 2003;65:489–500.View ArticleGoogle Scholar
- Stein AC, Uy JAC. Unidirectional introgression of a sexually selected trait across an avian hybrid zone: a role for female choice? Evolution. 2006;60:1476–85.PubMedGoogle Scholar
- Toews DPL, Hofmeister NR, Taylor SA. The evolution and genetics of carotenoid processing in animals. Trends Genet. 2017;33:171–82.PubMedView ArticleGoogle Scholar
- Blount JD, McGraw KJ. Signal functions of carotenoid coloration. In: Britton G, Liaaen-Jensen S, editors. Carotenoids. Basel: Birkhauser; 2008. p. 213–36.View ArticleGoogle Scholar
- Hey J. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics. 2004;167:747–60.PubMedPubMed CentralView ArticleGoogle Scholar
- Jankowski JE, Londoño GA, Robinson SK, Chappell MA. Exploring the role of physiology and biotic interactions in determining elevational ranges of tropical animals. Ecography. 2012;36:1–12.View ArticleGoogle Scholar
- Engler JO, Rödder D, Elle O, Hochkirch A, Secondi J. Species distribution models contribute to determine the effect of climate and interspecific interactions in moving hybrid zones. J Evol Biol. 2013;26:2487–96.PubMedView ArticleGoogle Scholar
- Chunco AJ. Hybridization in a warmer world. Ecol Evol. 2014;4:2019–31.PubMedPubMed CentralView ArticleGoogle Scholar
- Billerman SM, Murphy MA, Carling MD. Changing climate mediates sapsucker (Aves: Sphyrapicus) hybrid zone movement. Ecol Evol. 2016;6:7976–90.PubMedPubMed CentralView ArticleGoogle Scholar
- Taylor SA, White TA, Hochachka WM, Ferretti V, Curry RL, Lovette IJ. Climate-mediated movement of an avian hybrid zone. Curr Biol. 2014;24:671–6.PubMedView ArticleGoogle Scholar
- Parsons TJ, Olson SL, Braun MJ. Unidirectional spread of secondary sexual plumage traits across an avian hybrid zone. Science. 1993;260:1643–6.PubMedView ArticleGoogle Scholar
- Ruegg KC. Genetic, morphological, and ecological characterization of a hybrid zone that spans a migratory divide. Evolution. 2008;62:452–66.PubMedView ArticleGoogle Scholar
- Carling MD, Zuckerberg B. Spatio-temporal changes in the genetic structure of the Passerina bunting hybrid zone. Mol Ecol. 2011;20:1166–75.PubMedView ArticleGoogle Scholar
- Milá B, Toews DPL, Smith BT, Wayne RK. A cryptic contact zone between divergent mitochondrial DNA lineages in southwestern North America supports past introgressive hybridization in the yellow-rumped warbler complex (Aves: Dendroica coronata). Biol J Linn Soc. 2011;103:696–706.View ArticleGoogle Scholar
- Miller MJ, Lipshutz SE, Smith NG, Bermingham E. Genetic and phenotypic characterization of a hybrid zone between polyandrous Northern and Wattled Jacanas in Western Panama. BMC Evol Biol. 2014;14:227.PubMedPubMed CentralView ArticleGoogle Scholar
- Hudson RR, Coyne JA. Mathematical consequences of the genealogical species concept. Evolution. 2002;56:1557–65.PubMedView ArticleGoogle Scholar
- Engebretsen KN, Barrow LN, Rittmeyer EN, Brown JM, Moriarty Lemmon E. Quantifying the spatiotemporal dynamics in a chorus frog (Pseudacris) hybrid zone over 30 years. Ecol Evol. 2016;6:5013–31.PubMedPubMed CentralView ArticleGoogle Scholar
- Campagna L, Repenning M, Silveira LF, Fontana CS, Tubaro PL, Lovette IJ. Repeated divergent selection on pigmentation genes in a rapid finch radiation. Sci Adv. 2017;3:e1602404.PubMedPubMed CentralView ArticleGoogle Scholar
- Poelstra JW, Vijay N, Bossu CM, Lantz H, Ryll B, Müller I, et al. The genomic landscape underlying phenotypic integrity in the face of gene flow in crows. Science. 2014;344:1410–4.PubMedView ArticleGoogle Scholar
- Toews DPL, Taylor SA, Vallender R, Brelsford A, Butcher BG, Messer PW, et al. Plumage genes and little else distinguish the genomes of hybridizing warblers. Curr Biol. 2016;26:2313–8.PubMedView ArticleGoogle Scholar