Testing Hypotheses of Diversification in Panamanian Frogs and Freshwater Fishes Using Hierarchical Approximate Bayesian Computation with Model Averaging

Abstract

Most Neotropical frog and freshwater fish species sampled to date show phylogeographic breaks along the Pacific coast of the Isthmus of Panama, with lineages in Costa Rica and western Panama isolated from central Panama. We examine temporal patterns of diversification of taxa across this ‘western Panama isthmus’ (WPI) break to test hypotheses about the origin of species geographical distributions and genetic structuring in this region. We tested for synchronous diversification of four codistributed frog taxon-pairs and three fish taxon-pairs sharing the WPI break using hierarchical approximate Bayesian computation with model averaging based on mitochondrial DNA sequences. We also estimated lineage divergence times using full-Bayesian models. Several of our results supported synchronous divergences within the frog and freshwater fish assemblages; however, Bayes factor support was equivocal for or against synchronous or asynchronous diversification. Nevertheless, we infer that frog populations were likely isolated by one or multiple Pliocene–Pleistocene events more recently than predicted by previous models, while fish genetic diversity was structured by Pleistocene events. By integrating our results with external information from geology and elevational sea level modeling, we discuss the implications of our findings for understanding the biogeographical scenario of the diversification of Panamanian frogs and fishes. Consistent with the ‘Bermingham/Martin model’ (Molecular Ecology 1998, 7, 499–517), we conclude that the regional fish assemblage was fractured by processes shaping isthmian landscapes during the Pleistocene glaciations, including drainage basin isolation during lowered sea levels.

1. Introduction

By revealing the geographical histories of genetic lineages within species, phylogeography helps identify processes that influence the divergence, spread, and spatial-demographic fluctuations of lineages [1,2,3]. Single-species phylogeography surveys often reveal cryptic refugia and speciation events, for example [2,4,5]. By contrast, ‘comparative phylogeography’ presents a stronger approach for testing the generality of evolutionary patterns recovered within single species, often uncovering historical processes shaping species richness and genetic diversity [6,7]. Comparative phylogeography infers the histories of regional assemblages by testing for spatial and temporal phylogeographic congruence across codistributed species [1,8,9]. Congruent spatiotemporal divergences across taxa can indicate a shared history of responses to historical events [9,10].

Most comparative surveys in phylogeography are based on mitochondrial DNA (mtDNA) and are consequently subject to the caveats of using a single locus. Phylogeographic histories within individual species may be confounded by effects of stochastic variance in coalescent and mutational processes, natural selection, sex-biased dispersal, or reproductive success on genetic data [11,12,13,14,15]. Comparative studies guard against effects of such stochastic processes on inferences by searching for phylogeographic patterns replicated across multiple lineages [2,16]. This is analogous to improving evolutionary inferences through using multiple loci to estimate population parameters, multilocus models, or population/species trees across multiple gene genealogies [12,17,18]. Still, comparative approaches face major challenges to distinguish among potential scenarios that may have led to the observed patterns. Parameter estimates from mtDNA can suffer from low resolution, yielding wide confidence intervals that limit hypotheses testing [12]. The critical task of understanding the temporal framework of diversification to infer the historical biogeographical scenario underpinning spatial-genetic patterns can also be problematic, even if spatially congruent breaks are identified. Because variance in past effective population sizes (N_e), mutation rate, and generation time create high variation in gene tree topologies and their depths (correlated with timing of their divergences; e.g., [11,19]), what appear as multiple divergence events may actually reflect the obscured signal of a single event [20]. Thus, comparing independent gene-tree depths and divergence time estimates can lead to erroneous inferences of multiple vicariance events.

To guard against potential misinterpretations of biogeographical history, it is imperative to rigorously evaluate the timing of population divergences across species. Approximate Bayesian computation (ABC) methods for phylogeography allow accounting for potentially confounding stochastic coalescent effects while estimating parameters of phylogeographic datasets in such a way that sidesteps the need to calculate explicit likelihood functions. These “likelihood-free” methods permit simultaneously estimating demographic parameters and divergence times across populations by replacing the data with summary statistics, resulting in more efficient extraction of information, for example [21,22]. One broadly validated ABC method for comparative phylogeography, the MTML-msBayes (multi-taxa multi-locus msBayes) pipeline [22,23,24,25], permits comparing hypotheses in a parameter-rich hierarchical ABC (hABC) model using data from multiple codistributed population-pairs. This method has become widely used for testing explicit hypotheses about the timing of genetic divergences that have arisen during the assembly and diversification of regional species assemblages [26,27,28,29,30,31]. Based on simulations, some have cautioned that overly broad priors can result in low power to detect multiple divergences using MTML-msBayes [32]. However, this pitfall can be avoided by graphical sampling checks as well as comparing multiple MTML-msBayes model classes (priors) using ABC model averaging to account for uncertainty in model selection [24,33]). Other workers have also advocated leaving the summary statistics vector unsorted, as well as employing a Dirichlet-process prior as the hyperprior for the hyperparameter on the number of co-divergence events (Ψ) in the MTML-msBayes model and relying on this hyperparameter for inference [34,35,36]. However, our recent study using a cross-validation approach applied to simulated and empirical data, including the Panamanian datasets analyzed in the current paper, shows conclusively that choice of prior on Ψ has limited impact on inference, while re-sorting of summary statistics yields considerably more reliable results [37].

The Central American Isthmus has been an active area of phylogeographic research for nearly 20 years, with most studies focusing on inferring how, when, and why lineages colonized the isthmian landscape, diversified, and assembled into local communities [5]. A growing catalogue of species is being identified with common phylogeographic divergences across the Pacific-coastal lowlands of the Isthmus of Panama, with lineages in Costa Rica’s Golfo Dulce region and western Panama isolated from lineages in central Panama (Figure 1 and Figure 2; reviewed by [5]). This pattern has become known as the ‘western Panama isthmus’ (WPI) break [5]. Biogeographers interpret WPI breaks in frogs as resulting from multiple dispersals of different lineages into Central America across a continuous Miocene–Pliocene wet forest corridor, from both the north and the south [38,39,40,41], followed by vicariance caused by the development of tropical dry forest 12–4 million years ago (Ma; ‘dry forest hypothesis’; [42,43]). While frog fossils from the region are lacking, this scenario is supported by plant microfossils showing that modern dry forest species began appearing in central Panama up to mid-Pliocene [44,45]. Among freshwater fishes, the WPI break appears to reflect the effectiveness of rugged terrains at Soná and Azuero peninsulas (Figure 1) as biogeographical barriers to fish lineages, which arrived in multiple northwestward waves from South America since the Neogene. Bermingham and Martin [9] hypothesized that a corridor, formed by low sea level and tectonic uplift of the isthmus as the Cocos Ridge collided with the Panama microplate by ~7 Ma [46], permitted ‘early’ emplacement of freshwater fish species since Miocene–Pliocene. Subsequently, they argued that high eustatic sea levels of the mid-late Pliocene ~3.5–3 Ma [47] recreated a seaway over the isthmus, causing freshwater fishes to go extinct in central Panama ~3 Ma, before final isthmus emergence [9]. They hypothesized that the WPI break in fishes resulted from fragmentation due to drainage basin isolation during Pleistocene glacial periods over 2–0 Ma [9]. This ‘Bermingham/Martin model’ has also been invoked to explain phylogeographic patterns in pseudoscorpions, eels, catfishes, frogs, and other taxa, for example [5,48,49,50,51].

Clearly, previous studies generated or tested biogeographical models that provide convincing evidence for the WPI break. However, none of these studies explicitly accounted for genetic variance associated with coalescent processes, or among-lineage differences in demography or mutational rates, that may have influenced the variable patterns of genetic divergences of taxa associated with the WPI break ([9,42]; see Results). Thus, the temporal pattern and underlying mechanisms of diversification across the WPI break have not yet been fully explored. Additional studies are needed to evaluate whether the pattern of genetic structuring shared among frogs and freshwater fishes in the western Panama isthmus arose synchronously or asynchronously, and to determine the timing of historical events that have shaped community assembly and genetic divergence in this region.

The goal of this study is to examine the temporal patterns of diversification of frogs and freshwater fishes across the WPI break to test existing hypotheses about the origin of species distributions and genetic structuring in west-Pacific Panama. We use hABC simulations with ABC model averaging, and full-Bayesian divergence time estimation, to test the alternative hypotheses of a single, synchronous pulse of diversification, versus multiple pulses of diversification, in the frog and fish assemblages across this region. In doing so, we are able to directly test the hypotheses discussed above. Support for synchronous diversification correlated to Pliocene high seas ~3.5–3 Ma would agree with frog and fish WPI divergences arising from a common history of Neogene dispersal into western Panama followed by isolation wrought by fluctuating Plio–Pleistocene sea levels, and possibly other factors (e.g., dry forests isolating the frog populations). This would support the Bermingham/Martin model as a general model of diversification applicable to frogs and freshwater fishes. By contrast, an inference of one or multiple pulses of frog diversification <4 Ma would reject the dry forest hypothesis, whereas a pulse of fish vicariance events >2 Ma would reject the Bermingham/Martin model as an explanation for fish assemblage genetic structure across the western Isthmus of Panama.

2. Materials and Methods

2.1. Taxon Sampling and Molecular Data

We considered all phylogeographic studies of lower Central American (Costa Rica–Panama) taxa, including studies summarized in a recent comprehensive review [5], and subsequent publications. Taxa were included in this study if they met minimum criteria for our analyses: frogs or freshwater fish species/lineages codistributed across the study area; genetic divergence at the WPI break, with at least two individuals sampled on either side; and mtDNA gene sequences available from GenBank. We focused on the mtDNA locus because it presents a rapidly evolving and coalescing genomic region that often reflects population history [16], protein-coding mtDNA genes have been widely studied in Central American taxa permitting interspecific comparisons, and additional nuclear markers were not available across taxa (additional details in Discussion Section 4.2).

We obtained seven population-pair comparisons from eight taxa from our literature search (

Table 1. List of taxa used to examine temporal diversification patterns across the western Panama isthmus.

Taxon	mtDNA Genes	Length (bp)	ntotal (West/East)	π	Ti/Tv	% div	Source
Frogs
Agalychnis callidryas	16S, ND1	1149 (118, 1031)	28 (16/12)	0.0337	11.373	8.0	Robertson & Zamudio (2009)
Craugastor crassidigitus	cytb, cox1	1353 (714, 639)	13 (5/8)	0.0645	14.55	13.9	Crawford et al. (2007)
Dendropsophus ebraccatus	ND1, tRNAs	1877	18 (11/7)	0.0245	12.04	5.9	Robertson et al. (2009)
Engystomops pustulosus	cox1	564	15 (2/13)	0.0139	31.67	4.6	Weigt et al. (2005)
Freshwater Fishes
Andinoacara coeruleopunctatus	ATP6/8	842	10 (2/8)	0.0093	20.77	2.4	McCafferty et al. (2012)
Pimelodella chagresi	ATP6/8, cox1	1471 (842, 629)	14 (6/8)	0.0210	12.19	3.5	Bermingham & Martin (1998)
Roeboides occidentalis–R. guatemalensis	ATP6/8, cox1	1493 (842, 651)	11 (8/3)	0.0311	10.36	6.4	Bermingham & Martin (1998)

For multi-gene datasets, total alignment length is given in nucleotide base pairs (bp) followed by length of each gene in parentheses, in the same order genes are listed at left. The transition/transversion ratios (Ti/Tv; i.e., kappa parameters) are HKY + Γ + I (Hasegawa-Kishino-Yano [60] + gamma + proportion of invariant sites) model estimates, and percent genetic divergences (% div) are given as average pairwise HKY + Γ + I distances between population-pairs. π, nucleotide diversity [61].

). Taxa that met our criteria included the following four lineages of frogs: the red-eyed treefrog, Agalychnis callidryas [52]; the dirt frog, Craugastor crassidigitus [42]; the hourglass treefrog, Dendropsophus ebraccatus [53]; and the túngara frog, Engystomops (Physalaemus) pustulosus [51]. The following three freshwater fish lineages also met our criteria: the blue acara, Andinoacara coeruleopunctatus [54]; the Chagres catfish, Pimelodella chagresi [9,55]; and headstander tetras, Roeboides occidentalis–R. guatemalensis [9]. While these taxa are from different families and their broader geographical distributions differ to varying degrees, all species/lineages are codistributed in the study area and exhibit the WPI break (Figure 2), and most are thought to have expanded their ranges into Central America from southern source populations in eastern Panama or northern South America [9,39,40,56,57]. An exception to this is C. crassidigitus, which more likely colonized the region from an ancestral population in the north, that is, Nuclear Central America [41,58]. We excluded from our analyses two taxa inferred to exhibit the WPI break in previous studies, Brachyhypopomus occidentalis [9] and Synbranchus marmoratus [50]. In the case of B. occidentalis, available mtDNA sequences had insufficient sample sizes for estimating population parameters and coalescent modeling, which require ≥2 individuals from daughter lineages on either side of the break. Also, two studies conflicted as to whether this taxon showed the WPI break [9,59]. In S. marmoratus, the geographical location of the WPI break included a major north–south component of isolation across the Central Cordillera, thus the break was not as well confined to the Pacific coast as in other taxa we analyzed, making comparisons difficult.

A summary of sampling intensity and molecular characteristics of each dataset, and their published sources, is provided in Table 1. Our final database consisted of 109 sequences representing 10–28 individuals sampled from each of these seven species (2–16 individuals sampled from the west or east side of the break), permitting us to address the evolution of frog and freshwater fish species assemblages. Final mtDNA alignments for each species/lineage ranged from 564 to 1877 bp in length, and most (4/7) datasets contained cytochrome oxidase subunit 1 (cox1) ‘DNA barcode’ sequences. GenBank numbers and locality data for the samples used in each dataset are provided in Table S1. Sequence data were grouped a priori into populations east and west of the WPI break. We assumed Hasegawa-Kishino-Yano (HKY) [60] models with gamma-distributed rate variation (Γ) and an estimated proportion of invariant sites (I) wherever possible because (1) this model provides a more appropriate description of mtDNA sequence evolution (e.g., allowing recurrent mutations) than other nucleotide substitution models available in the MTML-msBayes pipeline (e.g., Jukes-Cantor or equal-input), (2) more complex models of sequence evolution are not implemented in MTML-msBayes, and (3) to make results comparable across analyses.

2.2. Genetic Diversity and Neutrality

We examined genetic diversity by calculating nucleotide diversity, π [61], for each dataset using the ‘obsSumStats.pl’ Perl script available in the latest version of the MTML-msBayes pipeline [25]. To further evaluate patterns of genetic divergence across the WPI break, we estimated average % divergence between clades for all population-pairs using distances corrected using the HKY + Γ + I substitution model in PAUP* version 4.0a [62]. In PAUP*, we specified Γ and I parameter estimates derived from Bayesian coalescent gene-tree analyses of each population-pair (below). Although selective neutrality is often a valid assumption for vertebrate organellar mtDNA, this was a key assumption of our analyses; thus, we tested for mtDNA selective neutrality using McDonald and Kreitman’s ([63]; MK) tests conducted in DnaSP version 5 [64]. Closely related, congeneric outgroup taxa were included in the MK tests for each population-pair based on phylogenies in the original studies, or large-scale phylogenetic studies (see Appendix A below for details).

2.3. Estimating Gene-Tree Depths Using Bayesian Dating

To independently assess variation in gene-tree depths across WPI break taxa, we estimated times to the most recent common ancestors (t_MRCAs) and their Bayesian credible intervals for each species/lineage using BEAST version 1.8.2 [65]. BEAST runs (100 million generations, sampled every 4000; ‘burn-in’ = 10%) started from UPGMA (unweighted pair group method with arithmetic mean) tree topologies, employed HKY + Γ + I site models, and used coalescent constant size tree priors [66]. We evaluated the fit of different molecular clock models and investigated clock-likeness of the data by comparing results of strict-clock and relaxed-clock (uncorrelated lognormal model; [67]) runs on each dataset. Due to uncertainty over frog mtDNA mutation rates, we specified a range of 0.1–1.3% lineage⁻¹ Myr⁻¹ (per million years) in the frog runs spanning rates of protein-coding mtDNA gene evolution documented in a broader mtDNA dataset from one of our focal taxa [51]. This ‘frog rate’ spans Macey et al.’s [68] Mongolian toad (Bufo bufo) rate of 0.69% lineage⁻¹ Myr⁻¹, which is commonly used to date patterns of herpetofaunal diversification [51]. Likewise, we specified a ‘fish rate’ of 0.17–1.4% lineage⁻¹ Myr⁻¹ in fish runs, a range spanning mutation rates estimated in 15 previous studies of teleost freshwater fish mtDNA (refs. in [69,70]). Rate ranges were supplied to the program as uniform priors. We summarized posterior distributions and ensured convergence and adequate effective sample sizes (all ESS >> 200) using Tracer version 1.5 (available at: http://beast.bio.ed.ac.uk/Tracer). In TreeAnnotator version 1.8.2, we summarized the posterior distribution of trees from each run by calculating a maximum clade credibility (MCC) tree annotated with median node ages from a sample of 5000 post-burn-in trees obtained using the ‘BEASTPostProc.sh’ script in PIrANHA version 0.1.4 [71].

2.4. Tests for Synchronous Diversification

We tested the hypotheses of a single, synchronous pulse of diversification versus multiple pulses of diversification of the seven population-pairs across the WPI break under a coalescent model, using the MTML-msBayes pipeline [25]. The parameters of the MTML-msBayes model are provided in

Table A1. Summary of generic MTML-msBayes model parameters for Y taxon/population pairs and their prior distributions (modified from [23,25]).

Parameter	Description	Prior Distribution
A. Hyperparameters
Ω	Dispersion index of divergence times, Var[τ]/E[τ] (ratio of variance to mean), across Y taxon/population-pairs	Set by Y, τmax, and Ψ
E[τ]	Mean multi-species (i.e., assembly-wide) divergence time across Y	Set by Y, τmax, and Ψ
Ψ	Number of distinct divergence times across Y	Discrete uniform [1, Y]
S	Vector of mutation rate scalars to accommodate among-locus mutation rate heterogeneity	Gamma (with shape and scale parameters α and 1/α, respectively)
α	Shape parameter of the gamma distribution	Uniform (1, 20)
B. Parameters
τ	Divergence time in units of 4N generations	Uniform (0.0, τmax)
θ = 4Nµ	Population size (mutation) parameter where N is average population size of two daughter populations, µ is the per-site mutation rate	Uniform (θmin, θmax)
θA	Ancestral effective population size (mutation) parameter	Uniform (0.01, θmax Nanc-max)
θA1, θA2	Ancestral effective population size (mutation) parameters for daughter populations 1 and 2 over ancestral period (from time τ to τB)	Uniform (0.01, θB1), (0.01, θB2)
θB1, θB2	Effective population size (mutation) parameters for daughter populations 1 and 2 over recent period (from time τB to 0.0, the present time)	Uniform (0.01 θ, 1.99θ), and where θB2 = 2θ − θB1
τB1, τB2	Times at which θA1 and θA2 experience exponential growth into populations with sizes of θB1 and θB2	Uniform (0.0, τ), with τB1 = τB2
M = Nm	Effective migration rate, where m is the probability of symmetric migration between pairs of daughter populations	Uniform (0.0, Mmax)
r	Per-locus mutation rate scalar	fixed
p	Per-locus ploidy and/or generation time scalar	p

. Because hyperprior ABC samplers like msBayes can inefficiently sample from the prior when overly broad prior bounds are used, biasing them towards inferring simultaneous diversification [32], we used ABC model averaging on candidate priors and visual checks for prior sampling efficiency to help overcome this issue [24]. Also, following recommendations for MTML-msBayes analyses of small numbers of taxon/population-pairs in an extensive simulation study [37], we ran MTML-msBayes with the default summary statistics and sorting algorithm, and we relied principally on the Ω hyperparameter (see below) for inference and hypotheses tests. Our MTML-msBayes analysis consisted of a four-step procedure similar to Hickerson et al. [24]. First, we developed four prior sets, or model classes (M₁–M₄), to compare using model averaging. Each model class consisted of one of two uniformly distributed priors on population divergence times (τ), ancestral population size (θ_A), and daughter population size parameters (θ_D;

Table 2. Prior model classes and results of tests for synchronous diversification using ABC model averaging in MTML-msBayes.

Prior	P(τ)	P(θD)	P(θA)	P(Mk|D)1000	Ω Mode	Ω 95% HPDs
WPI frogs (Y = 4)					0.0036	[0.000, 0.0565]
M1	~U(0, 1.75)	~U(0, 0.1)	~U(0, 0.25)	0.2666	–	–
M2	~U(0, 1.75)	~U(0, 0.1)	~U(0, 0.5)	0.4990	–	–
M3	~U(0, 1.75)	~U(0, 0.4)	~U(0, 0.25)	0.0000	–	–
M4	~U(0, 0.875)	~U(0, 0.1)	~U(0, 0.25)	0.2344	–	–
WPI fishes (Y = 3)					0.0017	[0.000, 0.0423]
M1	~U(0, 0.8)	~U(0, 0.1)	~U(0, 0.25)	0.3124	–	–
M2	~U(0, 0.8)	~U(0, 0.1)	~U(0, 0.5)	0.3766	–	–
M3	~U(0, 0.8)	~U(0, 0.4)	~U(0, 0.25)	0.0000	–	–
M4	~U(0, 0.4)	~U(0, 0.1)	~U(0, 0.25)	0.3110	–	–

Parameters are shown for four prior model classes (M_k) run for each of two analyses of Y population-pairs used to test for synchronous diversification across the western Panama isthmus (WPI) break. Prior models had varying τ, θ_D, and θ_A prior distributions, P(x) for random variable x, but assumed zero post-divergence migration. Approximate posterior probabilities, P(M_k|D)¹⁰⁰⁰, of each model are given based on 1000 accepted simulated draws from 20 million random draws from the four prior models, with that of the best-supported model underlined. Modal Ω hyperparameter estimates and their 95% highest posterior densities (HPDs) from model averaging over all four prior models are given in the first row of each section. Abbreviations: M, model class (candidate prior); P, probability (density); U, uniform distribution.

). Second, we obtained 5 million random (simulated) samples from each model class specified by a discrete uniform hyperprior distribution Pr(M_k) = 1/4. We visually checked for efficient prior sampling by conducting principal components analysis on 1000 prior draws from each model class in R. Third, we obtained the ABC joint posterior distribution using the default summary statistic vector (D) from MTML-msBayes and rejection sampling to identify the 1000 closest Euclidean distances between the observed summary statistics (D*) for the data and D_i calculated from 20 million random draws across all four priors (M₁–M₄). However, prior to rejection sampling, we rescaled the dispersion index of population divergence times Ω (=Var[τ]/E[τ]; the ratio of variance to the mean of the divergence times) and the mean assembly-wide divergence time, E[τ], from models with smaller upper θ_D prior bound values (frog and fish M₁, M₂, and M₄) to have the same coalescent units as the other model, M₃ [24]. Resulting estimates of Ω and E[τ] were weighted by Bayesian model averaging [24]. Last, we conducted hypothesis testing by comparing hyperposterior probability distributions of Ω estimates, to determine whether the data supported single or multiple diversification periods.

We conducted hypothesis testing by comparing posterior probabilities of the hyperparameter estimates under a ‘null’ scenario of asynchronous diversification (Ω > 0.01) versus the alternative of synchronous diversification (Ω ≤ 0.01) [22,23,28,30] and subsequently calculating B₁₀ Bayes factors under the parameter thresholds above while accounting for prior support for the hypotheses, using B₁₀ “weight of evidence” criteria in Jeffreys [72] and Kass and Raftery [73]. We estimated mean assemblage-wide divergence times by converting model-averaged E[τ] estimates (in coalescent units of 4N_ave generations) to absolute time (T_div) using the equation T_div = E[τ] × (θ_ave/μ), where μ is the assumed mutation rate per site per generation and θ_ave (per site) is the mean of the upper θ prior. Conversions used mutation rates equivalent to 0.7% lineage⁻¹ Myr⁻¹ and 0.785% lineage⁻¹ Myr⁻¹, the median rates of uniform ‘frog rate’ and ‘fish rate’ priors used in our BEAST analyses. We assumed an average generation time of 1 year, which was also used in previous studies of our focal taxa ([51,52,54]; Michael J. Ryan, pers. comm.), and we acknowledge that divergence time estimates are sensitive to generation times. Shell and R scripts used during our MTML-msBayes analyses are accessioned in Mendeley Data [74].

3. Results

3.1. Genetic Diversity and Neutrality

Clades examined here had considerably variable levels of genetic variation, with π values ranging 6-fold from 0.0093 in A. callidryas to 0.0645 in C. crassidigitus (Table 1); these same two taxa exhibited the lowest and highest levels of sequence divergence, respectively. There was also more variation in nucleotide diversity among frogs than among freshwater fishes. Likewise, genetic divergences based on HKY + Γ + I corrected distances between clades on either side of the break ranged from 4.6% up to a maximum of 13.9% among frog taxa, whereas that among fish species/lineages showed a more limited range from 2.4% to 6.4% (Table 1). Consistent with the assumptions of our analyses, MK tests did not reveal departures from selective neutrality within any of the mtDNA datasets (Fisher’s exact test p > 0.10; Appendix A).

3.2. Estimating Gene-Tree Depths Using Bayesian Dating

Estimated t_MRCAs for each of the seven population-pairs were consistent across multiple BEAST runs, which had ESS values >500 for nearly all parameters. Results also were similar across strict-clock and relaxed-clock runs; however, we only present the results of the strict-clock analyses because 95% highest posterior densities (HPDs) of ‘ucld.stdev’ (standard deviation of the relaxed clock) abutted zero, indicating that the data could not reject strict clock models. The mitochondrial MCC time trees had variable gene-tree depths (Figure 3), and geometric mean t_MRCA estimates (closer to peak likelihood values than the means) varied substantially across species/lineages, ranging from 1.71 Ma in A. coeruleopunctatus to 10.79 Ma in C. crassidigitus (Figure 4; frog geometric mean t_MRCA range: 3.78–10.80 Ma; fish geometric mean t_MRCA range: 1.71–4.88 Ma). Consistent with patterns of DNA polymorphism above, t_MRCAs of the fish lineages exhibited less variation with a much narrower region of overlap in their coalescence times (1.92–5.57 Ma) as compared with that of the frog lineages (3.88–14.39 Ma), although the frog t_MRCAs also overlapped substantially (Figure 4). Hereafter, geometric mean values are relied upon because they were closer to peak likelihood parameter estimates.

3.3. Tests for Synchronous Diversification

Projecting the observed data into principal components calculated from summary statistics randomly drawn from each model class showed that MTML-msBayes efficiently sampled our prior distributions (Figure S1). Observed data points also had similar positions across model classes, suggesting different priors were similarly good samplers for the data. Below we report results based on Ω posteriors derived from local linear regression [21].

We found that the modes of model-averaged posterior estimates of hyperparameter Ω were consistent with synchronous divergences within the frog and fish assemblages (frog Ω mode = 0.0036; fish Ω mode = 0.0017; Table 2; Figure 5). In agreement with this pattern, Ω = 0 fell within the credible intervals (95% HPD) of each of the model-averaged Ω distributions. The best-supported models also indicated synchronous divergence of the frog and the fish lineages based on posterior modal Ω values of 0.0013 and 0.0047, respectively. Despite this, model comparisons using Bayes factors for the ABC model-averaging results indicated that the data provided equivocal support for or against either model. Bayes factors were around a value of 1 for both the synchronous diversification model (frog B₁₀ = 1.14 for Ω < 0.01 versus Ω > 0.01; fish B₁₀ = 1.33 for Ω < 0.01 versus Ω > 0.01) as well as the asynchronous diversification model (frog B₁₀ = 0.88 for Ω > 0.01 versus Ω < 0.01; fish B₁₀ = 0.75 for Ω > 0.01 versus Ω < 0.01). Bayes factors were technically less than 1 for two or more divergences, but at best this provides very weak negative evidence for asynchronous divergence [72,73]. Likewise, posterior probabilities for the best-supported models were low (<0.5) and corresponding Bayes factors for synchronous diversification were weak, being approximately less than or equal to 1.

Based on model-averaging, the modes of assemblage E[τ] for the frog lineages and fish lineages were 0.106 (95% HPD = 0.056–0.139) and 0.054 (95% HPD = 0.022–0.079), respectively (Figure 5). When converted to absolute time using the median frog and fish mutation rates above, these peak-likelihood estimates correspond to modal dates of 3.09 Ma (95% HPD = 1.62–4.03 Ma) in the Pliocene to early Pleistocene for the frog divergences and 1.36 Ma (95% HPD = 0.56–2.021 Ma) in the early-mid Pleistocene for the fish divergences (Figure 5).

4. Discussion

4.1. Comparative Phylogeography of Panamanian Frog and Fish Assemblages

Determining whether species assemblages experienced shared evolutionary responses to historical geological and climate-change events is a central but challenging goal of comparative phylogeography [2,6,10]. Multi-taxon hABC methods for comparative phylogeography, for example, [22,23,24,25], permit testing explicit hypotheses about the timing of genetic divergences that have arisen during the assembly and diversification of regional species assemblages [26,27,28,29,30,31,33]. We used hABC coalescent models and full-Bayesian divergence dating to investigate temporal patterns of diversification of frogs and freshwater fishes across the WPI phylogeographic break (Figure 1 and Figure 2; [5]) to test two a priori hypotheses about the origin of species distributions and genetic structuring along the Pacific coast of Panama.

Several aspects of our results indicated that the frog and freshwater fish assemblages each experienced a synchronous pulse of diversification in this region during the late Neogene to Quaternary. This was supported by overlap in estimated t_MRCAs for the taxon-pairs (Figure 4), strong peaks in the Ω posteriors near zero (Figure 5) and credible intervals including zero (Table 2), and peak-likelihoods and credible intervals of assemblage E[τ]. By contrast, Bayes factor model selection indicated that the data provide only a marginal accumulation of evidence in favor of synchronous diversification and against the null scenario of asynchronous diversification. The msBayes approach has been shown to correctly reject synchronous diversification using only mtDNA and summary statistics such as those used herein (e.g., [23,24,25]), and visual checks suggested that our priors were efficient samplers of the data (Figure S1). Moreover, we avoided the problem of overly broad or narrow τ priors (e.g., in [26,32]), which can cause ABC samplers to explore parameter space exceeding saturation effects on mitochondrial genes [24], by matching the upper bounds of these priors to empirical estimates from the data. As a result, we conclude that additional genetic data from unlinked loci or additional species, or improved methods for hABC or Bayes factor estimation, are needed to more confidently assess the timing and number of events at the WPI break in these taxa. Nevertheless, the much narrower credible intervals of our assemblage divergence times relative to the Bayesian gene-tree depths inferred in BEAST (Figure 3 and Figure 4) indicate that accounting for coalescent processes and changes in population sizes through time in our models yielded much more precise estimated divergence times across the WPI break than were previously available. Assuming that peaks in Ω and E[τ] contain one or multiple clusters of population divergence events, and acknowledging limitations and caveats of our mtDNA data (see Introduction and Section 4.2 below), we use our Bayesian assemblage E[τ] estimates to draw broad conclusions about the timing of diversification in these two assemblages and conduct hypotheses tests. We also discuss implications of our results for understanding the historical biogeography and diversification of Panamanian frog and freshwater fish species assemblages.

The dry forest hypothesis predicts that the timing of diversification of frogs across the WPI break was marked by the transition from wet tropical forest in the middle Cenozoic to drier forest habitats along the Pacific coast of Panama by >4 Ma according to paleobotanical data [44,45], and up to 12 Ma based on mtDNA phylogenetic dating analyses [42,43]. Against these predictions, we inferred that frog diversification across the WPI break occurred more recently, around 3.09 Ma (Figure 3 and Figure 4). Credible intervals for model-averaged assemblage E[τ] from MTML-msBayes slightly overlapped the predicted interval, by only ~31 thousand years (kyr); thus, we do not reject the dry forest hypothesis outright. However, modal Ω and assemblage E[τ] estimates favor the conclusion that the fracturing of western and eastern Pacific Panama frog assemblage coincided with one or multiple historical events occurring after 4 Ma, during the late Pliocene–early Pleistocene. Differences between frog divergence times inferred by earlier workers and this study likely owe in part to differences in sampling and methods used. Crawford et al. [42] dated the isolation of Golfo Dulce Craugastor populations using analyses that combined interspecific samples and employed strict molecular-clock and nonparametric rate smoothing methods. Thus, unlike population divergence times we estimated in MTML-msBayes, earlier estimates of Golfo Dulce frog divergence times were based on gene divergences, which are more apt to overestimate the timing of diversification, for example [12].

While we cannot reject the dry forest hypothesis outright, our results are more consistent with the interpretation that Plio–Pleistocene sea level highstands limited frog dispersion and gene flow across the region, as predicted by the Bermingham/Martin model [9]. The Pliocene warm period was a time of reduced global ice volume, stable and warmer temperatures, and widely fluctuating sea levels 3.5–2.5 Ma [47,75]. Geological evidence from deep-sea sediments and passive continental margins shows that a Pliocene sea-level highstand +35 ± 18 m above present sea levels (a.s.l.) occurred 3.5–3 Ma [47,76] and swept over the emerging central Panama isthmus just prior to its final formation around ~2.8 Ma ([46,47,77]; but see [78]). Our divergence time estimates from MTML-msBayes place diversification in the frog assemblage around the time of this eustatic highstand (it falls within the Bayesian credible intervals completely) and match closest to its lower boundary (Figure 4).

To provide a post-hoc test of whether a Pliocene eustatic highstand could have isolated frogs across the WPI break based on external geographical data, we developed a geographic information systems (GIS) model of sea-level rise over the NASA Shuttle Radar Topographic Mission 250 m resolution digital elevation model (CGIAR-CSI; available at: http://srtm.csi.cgiar.org) using ArcMap version 10 (ESRI, Redlands, CA, USA). The model was constructed by tracing the 50 m, 100 m, and 250 m elevation contours at high resolution, similar to [5,79]. Assuming the upper limit of Cronin and Dowsett [47], the GIS sea-level model suggests that a highstand ~+50 m a.s.l. could have inundated the coastal plain between Golfo Dulce and Soná peninsula, the valley between Soná and Azuero peninsulas, and much of the east coast of Azuero peninsula, isolating it from central Panama (Figure 6). This model is conservative, given that lower Central American lands approached their modern subaerial extent by at least mid-late Pliocene if not earlier (see [78]), and the Pliocene Panama landscape was likely lower and more topographically homogeneous than the present-day landscape (reviewed by [5,44,80]). We hypothesize that the WPI break was subsequently reinforced by the mid-Pleistocene eustatic highstand ~550–390 ka, when sea levels were +22 m a.s.l. [77,81]. Relative to the dry forest hypothesis, this ‘marine vicariance hypothesis’ presents a testable and non-mutually exclusive hypothesis that explains how a eustatic highstand might have sparked multiple responses to a single period of sea-level rise and maintained any previous genetic divergence at the WPI break originally created by expansion of dry forests.

Within the Panamanian frog assemblage, lineage dispersal and persistence appear to be at least partly controlled by ecological requirements and degree of specialization, as well as life-history mode. For example, E. pustulosus is a generalist occurring in savannahs and open environments as well as humid and dry forests [40]. In our results, this species exhibited the youngest divergence across the WPI (Figure 3 and Figure 4), possibly reflecting greater capacity to disperse and establish populations in different habitats, or exchange genes with other populations. That said, the degree of genetic admixture among Panamanian populations of E. pustulosus is poorly understood, and has only previously been documented by nuclear allozymes [51]. Dendropsophus ebraccatus and A. callidryas are ecologically similar species ([40,53], refs. therein) and display broadly overlapping t_MRCAs slightly older than that of E. pustulosus but younger than that of C. crassidigitus (Figure 3). In turn, C. crassidigitus (like C. talamancae) is considered more ecologically specialized, preferring wet forest habitat more so than its dry-forest congeners (i.e., C. fitzingeri; [42]). That the C. crassidigitus gene tree extends farthest backward in time thus seems to suggest this species experienced long-term isolation and persistence in preferred habitats, rather than an ability to tolerate climatic or vegetational shifts. However, whereas all the other frog species are egg-layers, Craugastor dirt frogs are direct-developing species that readily reproduce, and their local population sizes can therefore be notoriously large [82]. Given the direct relationship between N_e and time to coalescence from coalescent theory, this contrast in life-history strategies would suggest that on one hand C. crassidigitus may have been superior at colonizing open niches or patches and spread throughout isthmian wet forest habitats more easily than the other taxa, while on the other hand its t_MRCA estimate may also be inflated due to large ancestral N_e [19].

Freshwater fishes typically exhibit higher levels of mtDNA sequence divergence and endemism in Atlantic relative to Pacific drainages of the Isthmus of Panama [9]. On the Pacific versant, the ~124 km-wide Gulf of Panama continental shelf (Figure 1) would have been broadly exposed as eustatic sea levels lowered 110–135 m during the Last Glacial Maximum (LGM, ~21 ka) and other glaciations, for example [75]. The above pattern likely owes to greater opportunities for dispersion and gene flow between drainages that coalesced over Pacific shelf areas during Pleistocene glacial periods [9,80]. Yet, this trend breaks down slightly in western Panama, where freshwater fishes show distinct associations between biogeography and phylogeography. Despite low endemism within fish provinces, geographical ranges of 14 fish species terminating at Soná and Azuero peninsulas (Figure 1) support a boundary between the Chiriquí and Santa Maria fish biogeographical provinces at the Río Santa Maria, Soná peninsula [80]. The WPI genetic break occurs in the same area in freshwater fishes (Figure 2C,D). Thus, historical processes importantly shaped Panamanian fish species turnover and diversification at the WPI break zone. Yet what mechanism best explains the isolation of fishes across the WPI break?

The Bermingham/Martin model predicts that WPI breaks in fishes have resulted from the isolation of fish populations west of Soná peninsula from those between Soná peninsula and El Valle volcano during glacial periods with sea-level lowstands since 2 Ma in the Pleistocene [9]. Consistent with this model, we found that the lower 95% HPDs for t_MRCA estimates for all three fish species/lineages fell within the predicted interval of diversification (Figure 3 and Figure 4). Peak posterior distributions from hABC model-averaging (e.g., Figure 5) also suggested that the focal fish species/lineages most likely diversified across the WPI break 1.36 Ma in the early Pleistocene, with Bayesian credible intervals ranging from early–late Pleistocene. This time period correlates best to the ‘Calabrian’ age (1.806–0.781 Ma; [83]), a time of 41-kyr periodicity of Pleistocene glaciations with drier and cooler-than-present conditions but less-extreme climatic oscillations than those following 800 ka [75,84]. Nevertheless, glacial periods vastly dominated the Calabrian to present, such that the Panama isthmus would have experienced many glacio-eustatic cycles but spent the majority of time since 1.8 Ma under glacial conditions with lowered sea levels and exposed continental shelf habitats. Eustatic sea-level curves give no convincing evidence that the oceans reached modern sea levels for any substantial period of time (e.g., >10–20 kyr) since the Calabrian, and the next eustatic sea-level highstand is not registered in the geological record until ~550–390 ka [76,81]. Geological patterns and processes are also consistent with decreased likelihood of fish dispersal across the WPI break zone since the Calabrian. In the break zone, the Pacific continental shelf becomes narrower (~0–40 km), tapering to the western Azuero peninsula coastline before being bisected by Cébaco and Coiba islands at the nearby Gulf of Montijo draining Soná peninsula (Figure 1). To evaluate the impact of lowered sea levels during Pleistocene glaciations on drainage connectivity in this area, we obtained a GIS model predicting paths of LGM paleo-drainages over modern bathymetry using ArcMap (courtesy of Peter J. Unmack, University of Canberra). The GIS model suggests that rivers draining to the west versus east of Soná peninsula did not anastomose over the continental shelf during the LGM (Figure 1), and possibly also preceding glaciations. Overall, our results combined with external environmental data suggest that a relatively stable geological setting at the Soná peninsula barrier has aided the historical isolation of drainage basins, maintaining fish lineage divergences at the WPI break during lower seas of the Calabrian to present.

While this study chiefly focuses on comparative phylogeographic modeling of codistributed Panamanian frogs and freshwater fishes, these assemblages also provide a system with ecological affinities (e.g., associations with freshwater pools and streams) and contrasts (terrestrial versus freshwater habits) making them suitable for assessing possible influences of the degree of ecological differentiation on phylogeographic structuring. One obvious ecological pattern among our results is that, at the deepest phylogenetic and ecological division among taxa in our study, frogs and fishes may have experienced largely independent pulses of diversification in the WPI break zone. Divergence events in these assemblages are not entirely statistically independent, given the overlapping credible intervals of their divergence times (Figure 3 and Figure 4), and our analyses treated these assemblages separately. However, the two-fold difference in the peak-likelihood E[τ] estimates suggest that frogs and freshwater fishes may have responded differently to the fluctuating paleogeographic, geologic, and climatic context of western Pacific Panama, for example [10]. We recommend that this hypothesis be tested using hABC models for comparative phylogeography improved by additional sampling (taxa and loci) and methodological improvements increasing the resolution of the models. However, under this scenario, dispersal of different lineages would be more strongly controlled by differences in ecological requirements among, rather than within, the two species assemblages. Thus, a potential biogeographical explanation for why these assemblages might have failed to register the impact of the same historical events is that the focal frog species/lineages, but not the fish lineages, colonized the Pacific coast of Costa Rica and Panama prior to the Pliocene eustatic highstand event, which left a greater signature on their genotypes than subsequent historical processes. This is supported by phylogenetic dating analyses of one frog lineage that we evaluated, the Craugastor fitzingeri group containing C. crassidigitus, which likely colonized the Central American Isthmus in the Eocene–early Miocene from a South American source population [58].

4.2. Caveats and Potential Limitations

4.2.1. Mitochondrial DNA

Our study makes several novel contributions, for example by providing the first hABC modeling study illuminating comparative phylogeography of both terrestrial and freshwater species assemblages from the Isthmus of Panama. There are limitations to an approach based solely on mtDNA, as also discussed in the Introduction. Therefore, we have provided a cautious interpretation of our results with the following limitations and caveats in mind.First, our hABC modeling analyses in MTML-msBayes accounted for coalescent stochasticity of the genetic data by allowing rates of lineage sorting to vary among population-pairs, but if multiple unlinked nuclear loci had been available to us, these may have allowed us to improve our accounting for the stochastic variance of coalescent processes and arrive at more accurate estimates of divergence times or other parameters, for example [12,85,86]. Even so, previous simulation tests of MTML-msBayes using pseudo-observed datasets from five taxon-pairs show that a single locus outperforms smaller multilocus datasets [25]. Significant improvements (i.e., lower root mean square error, RMSE) in estimating hyperparameters Ω and E[τ] only come when the total number of unlinked loci reaches ≥16, and additional increases in estimator accuracy require greater than 32 loci [25]. As a result, given that we analyzed mtDNA collected in studies conducted over a 14-year period (Table 1), it will likely take several more years before phylogeographic datasets from multiple Panamanian frog and freshwater fish taxa are sufficiently large to improve upon our results by scaling available datasets up to ≥16–32 loci. Moreover, if nuclear loci indicate substantial recombination, incomplete lineage sorting, or admixture across the WPI break, then these will pose additional analytical challenges and, particularly in the case of including migration parameters [22,25], are likely to require many more unlinked nuclear loci to develop models that perform well. One reason for this is that substantial migration can hinder the msBayes model from accurately inferring or rejecting simultaneous diversification, for example [22].

A second caveat to our results is that they may have been influenced by errors in the mutation rates supplied to the software programs. We conservatively identified ranges of plausible mutation rates for protein-coding mtDNA genes in frogs and fishes, and used these to set uniform priors on μ parameters in BEAST and MTML-msBayes. In each case, we allowed the software to estimate μ, through ABC or MCMC inference, and we allowed μ values to vary among population-pairs. Thus, our analyses accounted for potential differences in μ among species/lineages, especially when estimating gene-tree depths in BEAST. However, had more accurate estimates of mtDNA divergence rates been available, our estimates of species/lineage and assemblage divergence times could have been more accurate, and we could have been able to reduce the number of parameters in our models by fixing the mutation rates to accurate estimates for each species/lineage, mainly in the case of BEAST. If, however, mutation rates are broadly similar among taxa as our μ priors assume but the prior ranges were off by some particular factor (e.g., by half, or two-fold), then we would expect to have obtained similar results during our tests for simultaneous diversification under different μ priors, correcting by that factor. We also note that the original studies that produced the datasets we analyzed (see Table 1) estimated lineage divergence times using methods that infer gene divergence times from single loci. Their results are thus more likely to have overestimated the timing of lineage divergences across the WPI break than our methods, which only estimated population divergence times, for example [12].

A third caveat is that, while our datasets met the assumption of mtDNA neutrality (see Results and Appendix A.2.1 of Appendix A), we have made assumptions about species demographic histories in each of our analyses that are common to many studies, but that may have introduced error into the results. For example, in BEAST we set constant population size tree priors assuming the entire sample including both daughter populations for each population-pair evolved under the neutral coalescent with constant population size through time [19]. However, in MTML-msBayes, ancestral population size segments that evolved within each daughter population were assumed to experience exponential population growth to present. Additional demographic information and more flexible modeling approaches allowing different population-size transitions in different parts of the population trees for each population-pair, had they been available, may have permitted more accurate demographic models yielding more accurate parameter estimates. Nevertheless, although we allowed population size changes in our models, our main goal was not to infer historical demography of each species, at large.

4.2.2. Migration and “Secondary Contact”

A demographic model of “secondary contact” would include a period of pure genetic divergence, with two lineages initially diverging in the absence of gene flow and then experiencing relatively recent secondary contact. We did not include this kind of model because it is not available in the state-of-the-art comparative phylogeography software pipeline that we used. The MTML-msBayes implementation of the msBayes hABC pipeline that we used incorporates two models: a pure divergence model, or “isolation model”, and a migration model [25]. The latter model is an isolation-with-migration type model, in which migration parameters (Nm) are included for each of the Y taxon/species pairs in the analysis (Table A1; Figure S2), and migration is allowed to vary among each of these pairs from the time of initial divergence to the present (except this occurs backwards in time, since the software uses a coalescent model); thus, this model is similar but not equivalent to a scenario of secondary contact.

We also did not include migration parameters in our MTML-msBayes analyses. The reason for this is that there is essentially zero evidence supporting genetic patterns reflecting incomplete lineage sorting or migration across the western Panama Isthmus (WPI) break among datasets analyzed in our paper. In fact, only one study focusing on just one of the seven species/lineages that we evaluated previously inferred any admixture or migration between WPI lineages—in the case of the frog, Engystomops pustulosus, based on admixture signals in nuclear allozyme data [51]. In total, ~86% of the lineages that we evaluated in the present paper are completely sorted across the WPI break, with no mtDNA evidence for interpopulation gene flow across the break zone (Figure 2 and Figure 3). Under the MTML-msBayes migration model, all lineages diverge in the presence of gene flow, and their migration rates are drawn from the same uniform prior, for example parameterized with Nm ~U(0, 1.0) meaning that migration is allowed to vary up to one individual per generation [28]. Incorporating migration rates into multi-species analyses of co-diversification in the MTML-msBayes pipeline is thus inappropriate for our data, and would likely have had adverse effects on model performance. This is because using migration models would force the software to estimate additional migration parameters for species pairs for which there is no data signal that is informative for migration (e.g., alleles shared between diverging species/lineages), and this could theoretically result in incorrect estimation of population size (θ) and divergence time (τ) parameters.

The effects of migration on MTML-msBayes analyses are fairly well known, and our decision not to include migration in our models is supported by several previously published studies. In particular, in a paper presenting the updated MTML-msBayes version of the msBayes pipeline, Huang et al. [25] showed through simulation that hyperposterior results were less accurate when migration was included in the model (isolation with migration, not secondary contact), but the hyperposteriors, hence inferences, were improved when the correct migration model was specified. For example, simulating with migration when the data were known to be from a model with no migration (Nm = 0) yielded hyperposterior Ω (divergence time variance) and E[τ] (divergence time) estimates with much wider variance (RMSE), obscuring interpretation (see Figures 6 and 7 in [25]). Additionally, Bell et al.’s [28] recent comparative phylogeography study of codistributed frog species from the Australian Wet Tropics inferred patterns of population structure similar to those we see in frogs and fishes across the WPI break in Panama, with their Figure 1 showing genetic evidence for allele sharing or potential migration/admixture in only one (20%) of five frog species surveyed. They used Bayes factors [73] to compare MTML-msBayes models with and without migration, and their mtDNA-only results from analyzing data very similar to ours overwhelmingly supported isolation models [28]. Adding a small number of nuclear markers also did not cause migration models to fit better than isolation models for the majority of species in their dataset [28]. Given similarities between data and phylogeographical patterns in our study and Bell et al.’s [28], it seems reasonable to expect that comparing isolation and migration models using our data would yield similar results.

There are two additional points about migration of relevance to our MTML-msBayes analyses. First, the MTML-msBayes framework is not highly robust to migration, as migration can increase gene tree variance, such as variation in tree shape or depth [87], and in theory should reduce resolution of hABC inference. However, there is some evidence that MTML-msBayes is increasingly robust to errors in estimating the main hyperparameters of the model, Ω and E[t], with increasing amounts of data (numbers of loci), with or without migration. Still, simulations show that increases in estimator accuracy mainly occur when the number of independent, unlinked markers rises above 16 loci (see Figures 3A and 4A in [25]). There are no additional nuclear markers available for the codistributed taxa that we have studied in the present paper, and we are unaware as to whether any additional data will become available for these species in the foreseeable future. Thus, our paper is not only the first to provide a synthesis of the available data, but the data analyzed are also the best data available to address the questions and hypotheses that we focus on in our manuscript, using an appropriate framework for comparative phylogeography.

5. Conclusions

We conducted the first tests for simultaneous diversification of frogs and fishes at the WPI break in western Panama (reviewed by [5]) by analyzing mtDNA sequences using ABC methods accounting for uncertainty in model selection using modeling averaging [24,25]. Current findings show that, despite good prior selection and sampling properties (also supported for our datasets by simulations in [37]), hABC tests for synchronous diversification using mtDNA datasets remain challenging and difficulties can arise in conducting hypotheses tests based on Bayes factors when analyzing small numbers of taxon/population-pairs. However, our study also demonstrates that, in such cases, integrating assembly-wide divergence time estimates, E[τ], from hABC with external information from geology and elevation data can still generate novel biogeographical insights and hypotheses towards refining the existing biogeographical context for a region, even despite the limitations of mtDNA (see Introduction and Discussion Section 4.2.1). Recent advances in hABC analysis, including extensions to genome-wide single nucleotide polymorphism (SNP) data [88,89] and better inference through buffering multi-taxon divergence times [37], suggest an exciting period ahead for comparative studies of the diversification of vertebrate species assemblages in Panama and other Neotropical ‘hotspots’. We encourage future phylogeographic studies of the frog and fish taxa analyzed herein, and other Panamanian taxa, to build upon the present foundation by developing and analyzing datasets with increased taxon and character sampling, while capitalizing on these computational advances. A particularly fruitful way forward in future examinations of the WPI break would be to develop SNP datasets, for example using ddRAD-seq [90], and then test for co-demographic expansion and co-diversification using aggregate site frequency spectrum methods [89]. Such approaches should provide the increased power and resolution needed to more confidently identify the number and timing of divergence events that have shaped vertebrate diversification across the WPI break, and other multi-taxon genetic breaks along the Isthmus of Panama.