Understanding the evolutionary history of the species in a particular region provides insights into how that fauna was formed. Of particular interest to biogeographers is examining the impact a geographical barrier had in generating temporal genetic diversity among codistributed species. We examined the impact a major New World barrier, the Isthmus of Tehuantepec (IT) in southern Mexico, had on a regional bird fauna. Specifically, genetic data from 10 montane-forest bird taxa were analysed using approximate Bayesian computation (ABC) to test the hypothesis of simultaneous intraspecific diversification at the IT. Because effective population size (Ne) has the greatest impact on coalescent times, thereby affecting tests of divergence among codistributed taxa, we chose priors for both current and ancestral Ne using empirical estimates of theta. The ABC method detected two discrete diversification events. Subsequent analysis with the number of diversification events constrained to two suggests that four taxa diverged in an older event, with the remaining six diverging more recently. Application of a range of mutation rates from 2.0 to 5.0% Myr−1 places both events within the Pleistocene or Late Pliocene, suggesting that fluctuations in montane habitat induced by climate cycles and a late Pliocene seaway may have fractured this montane bird fauna. The results presented here suggest this avian fauna responded in a relatively concerted fashion over the last several million years.
A common yet elusive goal in historical biogeography and comparative phylogeography is determining whether codistributed taxa have experienced simultaneous diversification owing to a common barrier (Nelson & Platnick 1981; Avise et al. 1998; Bermingham & Moritz 1998; Avise 2000; Edwards & Beerli 2000; Arbogast & Kenagy 2001; Zink 2002). Estimates of the timing of vicariant events can be obtained by applying a molecular clock to species phylogenies. However, such estimates have large variances owing to differences in effective population sizes, mutational rates and mutational rate heterogeneity (Avise 1989; Hudson 1990; Arbogast et al. 2002; Rosenberg & Nordburg 2002). Because of this variation, temporally independent events can be difficult to distinguish from one another and may consequently be interpreted as a single event. Conversely, this variation can cause a simultaneous event to appear to have occurred at multiple times (Riddle & Hafner 2006). A method that can reduce or accommodate uncertainty owing to the stochastic nature of the evolutionary process is required if we are to understand the temporal nature of divergence across common barriers (Edwards & Beerli 2000; Knowles 2004).
Statistical phylogeographic methods offer a potential solution because they can account for the stochastic nature of the evolutionary process while allowing tests of simultaneous diversification. These methods can even incorporate biologically realistic scenarios such as different effective population sizes, expanding populations and various levels of gene flow, making them potent tools (Knowles & Maddison 2002; Knowles 2004; DeChaine & Martin 2006; Hickerson et al. 2006). In this study, an approximate Bayesian computation (ABC) method implemented in the recently developed statistical phylogeographic program MSBayes (Hickerson et al. 2006, 2007) was used to test for simultaneous intraspecific diversification of montane-forest birds at the Isthmus of Tehuantepec (IT) valley.
The IT (figure 1) in southern Mexico is a major geographical barrier to a variety of taxa (e.g. Parker et al. 1996; Peterson et al. 1999; García-Moreno et al. 2004). A valley approximately 250 m above the sea level and some 200 km wide at its narrowest point characterizes the IT (Barrier et al. 1998). The central valley is bounded by three mountain chains, the Sierra Madre Oriental and Sierra Madre del Sur to the northwest and the Chiapas-Guatemala Highlands to the southeast (figure 1). Montane organisms occurring on either side of the IT are currently isolated from one another by the intervening valley. Despite its current configuration, reshuffling of montane forests and their constituent faunas probably occurred multiple times during the last several million years (Wells 1983; Hooghiemstra et al. 2006), providing numerous opportunities for populations to diverge in isolation and for genetic admixture to occur when populations were reconnected. Thus, currently isolated taxon pairs might have originated at different times or perhaps during a single event. Molecular sequence data were used to test the hypothesis of simultaneous intraspecific diversification of 10 montane bird taxa having distributions that span the IT. The results provide an assessment of the temporal nature of intraspecific diversification in montane birds across this key geographical barrier.
2. Material and methods
Ten taxa with distributions spanning the IT were examined. Eight study species are Passeriformes, divided into six families: Cyanocitta stelleri (Corvidae), Certhia americana (Certhiidae), Piranga flava (Thraupidae), Myadestes occidentalis (Turdidae), Parula superciliosa, Basileuterus belli, Myioborus miniatus (Parulidae) and Atlapetes albinucha (Emberizidae). The two remaining taxa include a hummingbird, Lampornis amethystinus (Apodiformes: Trochilidae) and a woodpecker, Picoides villosus (Piciformes: Picidae). These 10 taxa are found primarily in pine–oak (Pinus and Quercus) and evergreen forests from 900 to 3500 m in elevation, and all are presumed to be resident taxa (Howell & Webb 1995).
(b) Molecular data
New and existing mitochondrial datasets were analysed (table 1: electronic supplementary material, appendix S1). New data were obtained as part of ongoing projects investigating the evolutionary history of montane-bird taxa in the region. Sequence data from the NADH dehydrogenase subunit 2 gene (ND2) were obtained from muscle tissues using standard molecular laboratory protocols for mitochondrial data. The primers L5215 (5′-TATCGGGCCCATACCCCGAATAT-3′: Hackett 1996) and H1064 (5′-CTTTGAAGGCCTTCGGTTTA-3′; Drovetski et al. 2004) were used in PCR and sequencing. Sequences were aligned and edited using Sequencher 4.7 (Gene Codes, Ann Arbor, MI, USA). Aligned sequences were checked for internal stop codons and indels to ensure against nuclear copies. Mitochondrial sequences from published sources were obtained from GenBank (table 1; electronic supplementary material, appendix S1). Datasets were grouped into populations occurring east and west of the IT.
(c) Estimates of the time-to-most-recent-common-ancestor
To provide an assessment of the variation in gene-tree depths, a Bayesian MCMC analysis was employed to estimate the time-to-most-recent-common-ancestor (TMRCA) and the confidence interval for each taxon. Analyses were conducted in BEAST v1.4.8 (Drummond et al. 2002; Drummond & Rambaut 2007). An HKY + G + I substitution model was assumed. Owing to the uncertainty of avian mtDNA substitution rates (Lovette 2004), we employed a range of 2.0–5.0% Myr−1 (Arbogast et al. 2002, 2006; Lovette 2004; Weir & Schluter 2008) and an uncorrelated exponential molecular clock (Drummond et al. 2006) to covert TMRCA to millions of years ago (Ma). A UPGMA topology was used as a starting tree, all other priors were set to default values. The analyses were run for 20 million generations with the first 25 per cent of the sampled points removed as burn-in. The mean TMRCA and its 95 per cent highest posterior density (HPD) were plotted from trees sampled every 1000 generations using Tracer 1.4 (Rambaut & Drummond 2007).
(d) Estimates of gene flow
Gene flow between populations on either side of the IT was estimated using IM (Hey & Nielsen 2004) for two taxa (M. occidentalis and P. superciliosa) with polyphyletic trees. Each analysis employed 10 heated chains and was run for a minimum of 3 × 107 generations with a burn-in of 5 × 105. Gene flow was set to be symmetrical between populations. Two independent runs were conducted on each dataset. Convergence was assumed when effective sample sizes (ESSs) were greater than 100 for each parameter (Kuhner & Smith 2007).
(e) Tests of simultaneous diversification
An ABC (Beaumont et al. 2002) method implemented in the program MSBayes (Hickerson et al. 2006, 2007) was employed to evaluate the relative strength of two hypotheses: simultaneous divergence versus non-simultaneous divergence. This program can simultaneously analyse multiple datasets while allowing each taxon and population to undergo independent demographic histories and has been shown to be effective over a range of conditions with a single locus (Hickerson et al. 2006, 2007). In addition, the method performs well with small sample sizes (e.g. one to three samples/population: Hickerson et al. 2007). The empirical range of Tajima's (1983) π was used to set the lower and upper limit of the uniform prior for the maximum (current) size of the population Θmax (where Θ = 4 * N * μ (μ is the per gene mutation rate)). The averaged value of Watterson's (1975) θ (Θw) and its ratio to the mean of π (Θw/π) was used to set the uniform prior for the maximum ancestral population (Θanc-max) to Θmax. The prior for number of possible maximum number of divergence events (¥) equalled the number of taxa (¥ = 10). Posterior distributions and their 95 per cent quantiles for number of divergence events (¥) were constructed using 1000 draws from 500 000 simulated replicates. Because the initial test suggested more than one event (i.e. mode of ¥ > 1.5 and Ω ≫ 0; see below) had occurred, a subsequent analysis was conducted to determine the number of taxa (with mean and 95% quantile) that diverged during each event and the timing (E(t)) of each of those events. This analysis was run with ¥ constrained to equal the number of divergent events detected in the unconstrained analysis, for example ¥ = 2. All other priors and settings were identical to the initial-unconstrained analyses. The mode of the mean divergence time across taxa (E(t)) was converted to years before the present (YBP) using the formula t(0.5Θmax/μ) and two mutation rates 2.0 and 5.0% Myr−1 (Arbogast et al. 2006). For these analyses, the transition to transversion ratio (table 1) was estimated from an equally weighted maximum parsimony tree using the likelihood option, and the HKY model, under tree scores in PAUP* v. 4.0b10 (Swofford 2002). Hickerson et al. (2007) recommends using both ¥ and Ω (=Var(τ)/E(τ)) to evaluate the relative strengths of each hypothesis, where a low value (approx. 0) of Ω and a ¥ = 1 suggest the data fit a simultaneous model; thus, we used both statistics to evaluate the relative support for each hypothesis.
(a) Molecular data
Sequence data from the ND2 gene (1041 bp) were obtained for eight taxa. GenBank numbers, museum numbers and localities are provided in electronic supplementary material, appendix S1. The L. amethystinus (Cortés-Rodríguez et al. 2008) dataset consisted of 354 base pairs for each individual. Samples sizes ranged from 2 to 55 individuals per population with an average of 14.7 (table 1). The number of variable characters and Ti/Tv ratio averaged 49 (14–121) and 11.4 (2.1–30), respectively (table 1).
ESSs were over 500 for nearly all statistics including TMRCA; thus, we assume the trees were well sampled. Mean values in YBP for the TMRCA ranged from 3.1 × 105 for P. superciliosa to 3.0 × 106 for M. miniatus. In general, there was a considerable variation and overlap in most but not all coalescent times for the 10 taxa (figure 2).
(c) Estimates of gene flow
Gene flow estimates (figure 3), in units of females/generation, between populations east and west of the IT were low for both M. occidentalis (0.05) and P. superciliosa (0.18); therefore, we assume the lack of unsorted trees in these two taxa is the result of incomplete lineage sorting and not due to ongoing gene flow (Palsbøll et al. 2004). Independent runs converged on similar values for all parameters and ESS values for all parameters were greater than 1000; thus convergence upon stationarity was assumed (Hey 2005; Kuhner & Smith 2007).
(d) Test of simultaneous diversification
The ABC method detected two divergence events among the 10 taxa. The upper and lower limits for the uniform prior of Θmax were set using the empirical range of π (2.9–33.8). Because the ratio of the means of Θw and π was close to one (13.2/14.2 = 0.93), we assumed the maximum-ancestral population was equal in size to Θmax; thus Θanc-max = Θmax. The initial ABC analysis (figure 4) detected non-simultaneous divergence (¥ = 1.8: 1–4 and Ω = 0.23) among the 10 taxa. The analysis with number of divergence events constrained to two (¥ =2) suggested six (three to nine) taxa diverged in the more recent of the two events, with the remaining four (one to seven) diverging in the older one. As currently implemented, MSBayes will not sort taxa into each of the divergence events; instead we used relative depths of TMRCA (figure 2) to rank taxa. Using TMRCA suggests P. flava, L. amethystinus, M. miniatus and C. stelleri or C. americana split during the older event, with the remaining taxa diverging in the more recent one. The mode divergence times (E(t)) for the two events were 0.08 (0–0.21) and 1 (0.42–2.02). When converted to YBP, these values become 130 000 (0–341 000) and 1 623 000 (682 000–3 247 000) YBP using 2.0% Myr−1 and 52 000 (0–136 000) and 649 000 (273 000–1 312 000) using 5.0% Myr−1 (figure 5).
Comparative biogeographic studies have usually assumed that differences in coalescence times signal differences in divergence times, thus rejecting simultaneous diversification (Avise 1992; Arbogast & Kenagy 2001). However, coalescent theory predicts a large variation in the depths of gene trees even when multiple codistributed taxa have responded to the same barrier (Hudson 1990). Inspection of TMRCA of the 10 taxa reveals large (as predicted by theory), and mostly overlapping, variation in coalescent times, but with obvious gaps. For example, the 95 per cent HPD of M. miniatus and P. flava do not overlap with that of P. superciliosa and M. occidentalis. Thus, variation in TMRCA suggests the possibility of more than one divergence event. The question then becomes whether this variation is the result of a single or multiple divergence events. Examination of TMRCA alone does not allow the determination of the number of divergence events. We employed a recently developed phylogeographic statistical tool to determine the number of diversification events that occurred among the 10 taxa.
The ABC method used in this study tests for simultaneous diversification among codistributed taxa, while at the same time allowing the population histories (e.g. effective population sizes, population growth, etc.) to vary across taxa. Because effective population size (Ne) has the greatest impact on coalescent times and consequently on tests of divergence times among codistributed taxa, it was important that the prior for Ne approximates empirical values. We chose prior values for current and ancestral Ne using values obtained from the data. In this way, the Ne priors encapsulate the range of observed values. The ABC method detected two pulses of divergence among the 10 montane bird taxa, a result consistent with the variation in TMRCA. Application of a range of mutation rates places the timing of these two events at 130 000 (0–341 000) and 1 623 000 (682 000–3 247 000) YBP using 2.0% Myr−1 and 52 000 (0.0–136 000) and 649 000 (273 000–1 312 000) using 5.0% Myr−1. It is important to note that these dates are the population division and will be more recent than TMRCA estimates. Thus, both events, including their confidence intervals, fall within the Pleistocene and Late Pliocene.
Climate change during the Pleistocene and Late Pliocene had a profound impact on species (Hewitt 2000). Montane forests underwent repeated altitudinal migrations up- and down-slope during this period (Wells 1983; Hooghiemstra et al. 2006). It is likely that there existed multiple successive opportunities for populations to diverge in isolation. The amplitude of each climate oscillation, at least for recent cycles, appears to be of similar magnitude (Petit et al. 1999; therefore, when forests were reconnected during each successive glacial cycle, evidence of a previous divergence might have been erased (Zink et al. 2004). Owing to this fact, we might expect only the last glacial event to leave its genetic signature upon codistributed taxa. These data suggest that an Early Pleistocene or Late Pliocene divergence maintained its evolutionary independence over subsequent glacial cycles.
It is common in biogeographic studies to assume that climate-driven habitat shifts are responsible for shaping genetic patterns when events date to the last two million years (e.g. Weir & Schluter 2004). Although population division owing to habitat shifts is plausible in this system, it is inferred from a casual temporal relationship (as in most studies), and alternatives should be considered. For example, a marine transgression may have crossed the isthmus valley at the Plio–Pleistocene boundary (Campbell 1999; Mulchay et al. 2006; but see Miller et al. 2005). The postulated timing (approx. 2 Ma) of this seaway falls within the confidence interval of the deepest divergence event (figure 5). Thus, multiple mechanisms (habitat shifts, seaways), operating at different times, may be responsible for fracturing this montane bird community. Vicariance owing to a seaway, rather than habitat shifts, could explain why the oldest divergence event maintained its evolutionary independence over later glacial cycles. Finally, different natural history traits (diet, canopy versus terrestrial foragers, etc.) may play a role in generating two events. However, uncertainty in the exact number of taxa that were split by each event precludes rigorous analysis.
Several phylogeographic studies of montane taxa exist in the region, providing an opportunity to examine general patterns and processes across taxa. For example, Sullivan et al. (2000) conducted a comparative study of two highland rodents spanning the isthmus region and detected phylogeographic patterns indicative of an isthmus vicariant event in both taxa. Overlapping variation in sequence divergence suggests the possibility that both taxa split simultaneously in the Pleistocene or Late Pliocene (Sullivan et al. 2000). Estimates of divergence times in two codistributed snakes suggest a Mid-Pliocene event, although confidence intervals overlapped with the earliest event detected in this study (Castoe et al. 2009). Finally, additional bird taxa (not included here because they employed different gene regions or occur in different habitat) show phylogeographic breaks at the IT. For example, the common bush-tanager (Chlorospingus ophthalmicus) and the emerald toucanet (Aulacorhynchus prasinus) show recent phylogeographic breaks (García-Moreno et al. 2004; Puebla-Olivares et al. 2008) that are roughly consistent with the recent divergence event detected in this study. A deeper and perhaps older genetic division was recovered between populations east and west of the isthmus in the Buarremon brunneinucha complex (Cadena et al. 2007; Navarro-Sigüenza et al. 2008).
Understanding the temporal nature of diversification across a shared barrier is a fundamental yet challenging goal in evolutionary analysis. Documenting the evolutionary divisions within extant faunas will provide insights into how communities have responded to, and been shaped by, historical events (Brooks & McLennan 1991; Zink 2002; Cavender-Bares & Wilczek 2003). ABC suggests that montane bird fauna occurring on both sides of the IT was fragmented during two discrete pulses of divergence. These two events occurred during the Pleistocene or Late Pliocene, implying climate change and possibly a marine transgression were responsible for shaping intraspecific diversity in the region of the isthmus. Finally, we recognize the limitations of using a single locus to test hypotheses of divergence as well as date those events (Edwards & Beerli 2000; Carstens & Knowles 2007). However, these results provide the first comparative genetic-based assessment of the temporal nature of diversification in this montane bird community and a working hypothesis for further study in the IT region.
We thank A. Jones, R. M. Zink, G. Barrowclough, J. Johnson, P. Unmack, J. Bagley and J. Lai for valuable comments. P. Sweet, B. Smith, J. Lai and G. Groth provided assistance with tissues and laboratory work. M. Hickerson assisted with the use of MSBayes. P. Escalante facilitated fieldwork in Mexico. M. Mika, B. Smith, M. Gurrola, J. Lai and R. Bargiel provided assistance with field collecting. Finally, we thank the various academic institutions (electronic supplementary material, appendix S1) and people that provided tissues. Fieldwork was supported by the Frank M. Chapman memorial fund. This work is a contribution from the Monell Molecular Laboratory and the Cullman Research Facility in the Department of Ornithology, American Museum of Natural History and has received support from the Lewis B. and Dorothy Cullman Program for Molecular Systematics Studies, a joint initiative of the New York Botanical Garden and the American Museum of Natural History, the Sackler Institute of Comparative Genomics and the L.J. and L.C. Sanford funds. Analyses using IM and BEAST employed the Computational Biology Service Unit from Cornell University, which is partially funded by Microsoft Corporation.
- Received February 17, 2010.
- Accepted March 31, 2010.
- © 2010 The Royal Society