Coordinated Dispersal and Pre-Isthmian Assembly of the Central American Ichthyofauna.

We document patterns of coordinated dispersal over evolutionary time frames in heroine cichlids and poeciliine live-bearers, the two most species-rich clades of freshwater fishes in the Caribbean basin. Observed dispersal rate (DO) values were estimated from time-calibrated molecular phylogenies in Lagrange+, a modified version of the ML-based parametric biogeographic program Lagrange. DO is measured in units of "wallaces" (wa) as the number of biogeographic range-expansion events per million years. DO estimates were generated on a dynamic paleogeographic landscape of five areas over three time intervals from Upper Cretaceous to Recent. Expected dispersal rate (DE) values were generated from alternative paleogeographic models, with dispersal rates proportional to target area and source-river discharge volume, and inversely proportional to paleogeographic distance. Correlations between DO and DE were used to assess the relative contributions of these three biogeographic parameters. DO estimates imply a persistent dispersal corridor across the Eastern (Antillean) margin of the Caribbean plate, under the influence of prevailing and perennial riverine discharge vectors such as the Proto-Orinoco-Amazon river. Ancestral area estimation places the earliest colonizations of the Greater Antilles and Central America during the Paleocene-Eocene (ca. 58-45 Ma), potentially during the existence of an incomplete Paleogene Arc (∼59 Ma) or Lesser Antilles Arc (∼45 Ma), but predating the GAARlandia land bridge (∼34-33 Ma). Paleogeographic distance is the single best predictor of DO. The Western (Central American) plate margin did not serve as a dispersal corridor until the Late Neogene (12-0 Ma), and contributed relatively little to the formation of modern distributions.

The powerful central insight of historical biogeography is that Earth history events impose common influences on the evolutionary diversification of multiple species, higher taxa, and portions of whole biotas (Rosen 1978; Humphries and Parenti 1986, 1999; Wiley 1988; Crisci et al. 2003). Under the vicariance paradigm, geographic distributions are explained as arising from a common history of geographic-range fragmentation, and dispersal is treated as an idiosyncratic and uncoordinated process producing no clear pattern (Zink et al. 2000; Donoghue 2008, 2011).

A compliment to the vicariance paradigm explains concordant geographic ranges as arising from the erosion of barriers to dispersal and gene flow. This process has been referred to as “geodispersal” (Lieberman and Eldredge 1996; Lieberman 1997; Upchurch 2008), or “barrier removal” (Johnstone 1986; Júnior et al. 2009). Geologically coordinated patterns of dispersal involving multiple taxa have been implicated in the formation of diverse paleontological and neontological assemblages, including Paleozoic trilobites in Gondwana (Meert and Lieberman 2004; Couvreur et al. 2011) and Laurasia (Stigall Rode and Lieberman 2005), Triassic ptycholepid fishes in Pangaea (Mutter 2011), amphisbaenian reptiles in Western Gondwana (Hembree 2006), pantropically distributed annonacean plants (Couvreur et al. 2011), freshwater fishes of tropical South America (Albert and Crampton 2010; Albert and Carvalho 2011), and many continental taxa during the Great American Biotic Interchange (Webb 1976, 1991).

Coordinated patterns of dispersal in multiple taxa are also expected from common exposure to geologically persistent dispersal vectors (De Queiroz 2005; McGlone 2005; Cowie and Holland 2006; Samonds et al. 2012). For example, stable, long-term patterns of atmospheric and oceanic circulation (i.e., prevailing winds and sea currents), coupled with the position of geological terranes in time and space, have been implicated in the formation of dispersal-assembled faunas of the Greater Antilles (Matthew 1919; Darlington 1938; Hedges et al. 1992; Censky et al. 1998; Iturralde-Vinent and MacPhee 1999; Glor et al. 2005; Iturralde-Vinent 2006; Heinicke et al. 2007; Michelangeli et al. 2008; Pyron 2014; Lewis et al. 2015), and Central America (McCafferty 1998; Smith and Bermingham 2005).

The Greater Antilles and Central America (henceforth GACA) region has figured prominently in studies of how vicariance, geodispersal, and coordinated dispersal affect the formation of insular biotas (Rosen 1975, 1978; Hedges et al. 1992; Losos and Schluter 2000; Trejo-Torres and Ackerman 2001; Graham 2003; Losos and Ricklefs 2009). The geological history of the Caribbean plate is incompletely understood and the timing of many important tectonic events is poorly constrained. However, at least three major geological events (and the geographic structures they produced) have been hypothesized to affect the distributions of species among landmasses of the Caribbean basin. These events are the formation of: 1) the Upper Cretaceous to Paleocene (94–63 Ma) Caribbean Large Igneous Province and associated movement of the Proto-Antilles arc between South America and Central America (Holden and Dietz 1972; Malfait and Dinkelman 1972; Mann 1999; Pindell and Kennan 2009; Loewen et al. 2013), 2) the Eocene–Oligocene (ca. 34–33 Ma) Greater Antilles-Aves Ridge (GAARlandia) land bridge between northern South America and the Greater Antilles (Iturralde-Vinent and MacPhee 1999; Escalona and Mann 2011), and 3) the Middle Miocene to Pleistocene (12–3 Ma) uplift of the Isthmus of Panama (Bacon et al. 2015; Montes et al. 2015).

Freshwater fishes have long been used as empirical models in the study of GACA biogeography, as they exhibit limited ecophysiological tolerances to salty marine waters and low capacities to traverse marine barriers (Myers 1949, 1966; see Sparks and Smith 2005). The GACA ichthyofauna is composed of about 600 species distributed across the Eastern (Antillean) and Western (Central American) margins of the Caribbean plate. The majority of these species are secondary freshwater fishes, especially Cichlidae and Cyprinodontiformes, which exhibit some tolerance to brackish waters, and which can occasionally cross narrow marine barriers (Myers 1949, 1966). The GACA ichthyofauna is dominated by clades with Neotropical distributions, including 454 species (86%) in 20 clades with South American origins (Matamoros et al. 2012, 2015). The taxonomic composition of the GACA ichthyofauna is a filtered version of the Neotropical fish fauna as a whole, including about one-third of the orders and families, and about 10% of the species (Chakrabarty and Albert 2011). The two most species rich clades in the GACA region are secondary freshwater fishes of Neotropical origins, heroine cichlids (132 species) and poeciliine live-bearers (204 species), which together comprise about 74% of the regional species pool (Matamoros et al. 2012, 2015).

Parametric biogeography is an emerging field at the interface of macroevolution and historical biogeography that models the history of geographic range evolution within explicitly phylogenetic and spatial contexts (Ree and Sanmartin 2009; Sanmartín 2010). Current methods to estimate macroevolutionary dispersal and extinction rates use maximum-likelihood or Bayesian methods and a dispersal–extinction–cladogenesis (DEC) model of species range evolution (Ree et al. 2005; Ree and Smith 2008; Ree and Sanmartin 2009; Smith 2009; Buerki et al. 2011; Goldberg et al. 2011; Ronquist and Sanmartín 2011; Landis et al. 2013; Wood et al. 2013). A shortcoming of existing implementations of the DEC model, however, is the requirement for the user to input dispersal rates (Buerki et al. 2011; Ronquist and Sanmartín 2011; Landis et al. 2013; Matzke 2013a, 2013b). A preferable approach is to use data from time-calibrated phylogenies, and the geographic distributions of extant species, to estimate empirical rates of macroevolutionary dispersal.

Here we investigate alternative hypotheses for the formation of GACA freshwater fish assemblages using parametric biogeographic methods to estimate historical rates of coordinated dispersal in heroine cichlids and poeciliine live-bearers. We develop and implement new methods for: 1) measuring macroevolutionary dispersal rates in the context of a landscape with changing paleogeographic interconnections; 2) estimating expected rates of dispersal from paleogeographic considerations of area, river discharge volume, and distance; and 3) comparing observed and expected dispersal rate values to identify the relative contributions of these paleogeographic agencies.

Materials and Methods

Phylogenetic Methods and Divergence Time Estimates

We acquired sequence data from published sources available in GenBank. For cichlids these data include three nuclear (RAG1, 1475 bp, RAG2, 850 bp, and S7 intron 1550 bp) and four mitochondrial (16S, 550 bp; CytB, 1100 bp; ND4, 647 bp; and COI 542 bp) genes, for a total of 4662 bp in 164 cichlid species, including 97 heroine species, and 8 percomorph outgroups (Martin and Bermingham 1998; Farias et al. 2001; Concheiro-Pérez et a. 2007; Říčan et al. 2008; Smith et al. 2008; Lopez-Fernandez et al. 2010; Chakrabarty and Albert 2011). For poeciliines these data include two nuclear (RAG2, 954 bp and S7 intron, 712 bp) and four mitochondrial (16S, 596 bp; CytB, 1137 bp; ND4, 673 bp; and COI, 590 bp) genes, for a total of 3778 bp in 103 species; Schartll 1994; Lydeard et al. 1995; Hrbek et al. 2006). All gene sequences were independently aligned using MAFFT v. 5.3 under default parameters (Katoh et al. 2002, 2005). Optimum partitioning schemes and substitution models were estimated in PartitionFinder (Lanfear et al. 2012). Analyses started with a fully partitioned data set (by gene and by codon position in protein-coding sequences) and the best-fit partitioning scheme and nucleotide substitution models were selected based on corrected Akaike Information Criterion (AICc; Akaike 1974).

We performed Bayesian analyses in B v.1.47 (Drummond et al. 2006; Drummond and Rambaut 2007) to simultaneously estimate tree topology and lineage divergence times using a lognormal relaxed molecular clock, and prior paleontological constraints. B analyses were composed of two independent runs, each starting from a random tree, using birth–death priors for branching times, and sampling parameter values and single trees every generations. The two runs were combined using L v1.7.5. Model parameter values were inspected for stationary convergence in Tracer v1.5. To ensure adequate mixing of the Markov Chain Monte Carlo, we inspected parameter estimates for effective sample sizes. The initial generations were discarded as burn-in. Posterior probability densities were summarized as a 95% maximum clade credibility tree using T v1.7.5. Posterior probabilities at nodes, means, and 95% highest posterior density estimates of divergence times were visualized in FigTree v1.4.0.

We constrained absolute divergence times by using lognormal calibration priors based on dated fossils. The dates were assigned using geological records from the oldest known fossils for clades of cichlids and poeciliids. Cichlid fossils include: †Proterocara argentina (mean: 2.0, SD: 1.0, offset: 49.0, calibration: Crenicichla + Teleocichla) and †Gymnogeophagus eocenicus (mean: 2.0 SD: 1.0 offset: 49.0 calibration: Gymnogeophagus) from the Paleocene–Eocene Lumbrera and Maiz Gordo Formations in northern Argentina (Malabarba et al. 2006, 2010); †Tremembichthys garciae (mean: 2.0, SD: 1.0, offset: 33.9 calibration: Cichlasomatinae) from the Oligocene Tremembé Formation in southern Brazil (Malabarba and Malabarba 2008); †Aequidens saltensis (mean: 2.5, SD: 1.0, offset: 12.0, calibration: Aequidens) from the Miocene Anta Formation in Argentina (Bardack 1961); and †Nandopsis woodringi (mean: 2.5, SD: 1.0, offset: 6.0, calibration: Nandopsis) from the Miocene Lascahobas Formation in Haiti (Cockerell 1923). Offset values for the minimum age of cichlids were set at 95 Ma (mean: 2.0, SD: 1.0), as the approximate age of the oldest acanthomorph fossil taxon †Polymixia sp. (Patterson 1993). Poeciliid fossils include: 1) an undescribed species (mean: 2.0, SD: 1.0, offset: 40.0 calibration point: clade Poeciliidae) from the Salta Formations of northwestern Argentina estimated at about 40 Ma (Pascual et al. 1981) and †Prolebias (mean: 2.0, SD: 1.0, offset: 30.0, calibration point: clade Aphanius) from the middle Oligocene of Europe dated to about 30 Ma (Hrbek and Meyer 2003). Mean, SD, and offset represent a set of parameters for the lognormal prior resulting in a distribution curve with lower and upper bounder of absolute ages.

Parametric Biogeography

Our goal was to estimate observed (i.e., empirical) dispersal rates ( directly from time-calibrated molecular phylogenies on a geologically informed paleogeographic landscape. To accomplish this goal we modified the Python code of the parametric biogeography program L (Ree and Smith 2008; Ree and Sanmartin 2009), to search for multiple independent dispersal rate estimates. This modified version of L is referred to here as L + (available online at https://code.google.com/p/lagrangeplus/). L + differs from the original L program in treating each dispersal rate parameter in the dispersal matrix as a separate value to be estimated. The dispersal scalar is removed from the process, and the dispersal matrix is then interpreted as a set of absolute (not relative) dispersal rates. The models of area connectivity through time (Fig. 1, center panel) employ eight independent dispersal rate parameters among five geographic areas and three time intervals.

Figure 1.

Estimating observed dispersal rates ( across the Caribbean Plate margins in L+. Paleogeographic reconstructions (left) and models of area connectivity through time (center), represented as a dispersal rate matrices (right). Dispersal rate parameters (Greek letters) are estimated in a ML framework employing a Dispersal–Extinction–Cladogenesis (DEC) model of geographic range evolution and empirical time-calibrated molecular phylogenies. Time intervals (I–III) encompass major geological events thought to have affected vicariance and dispersal across the Caribbean plate (see text). Paleogeographic maps and data from Pindell and Kennan (2009) and Blakey (2011). Areas: A, Cis-Andean; B, Trans-Andean; C, Panama; D, Nuclear Central America; and E, Greater Antilles.

Sensitivity of the L+ ancestral area estimates to model assumptions was assessed using several alternative analytical methods implemented in B (Matzke 2013a, 2013b, 2014). These analyses differed in the model of biogeographic range evolution (i.e., DEC and DECj), and in how parameters were used to calculate ancestral areas using the L, D-L, B-A-L analytical methods. The D-L model uses likelihood not parsimony to estimate ancestral areas, and is therefore referred to as D-L. The B model is a likelihood interpretation of the Bayesian model implemented in the B program of Landis et al. (2013).

Dynamic Paleogeographic Model of the Caribbean Plate

The Caribbean plate occupies approximately 3.2 million km2 between the North American, South American, Nazca, and Cocos plates. The Eastern margin of the Caribbean plate is marked by the Antillean Arc with an emergent (sub-aerial) land area of about 220,000 km2, taken as the sum of the areas of modern Cuba, Jamaica, Puerto Rico, Hispaniola, Cayman Islands and Lesser Antilles. The Western margin of the Caribbean plate is marked by Central America, extending between the Isthmus of Tehuantepec in southern Mexico (∼19° North latitude) and the Isthmus of Panama (∼7° North latitude), with an emergent land area of about 2,368,000 km2 (Chakrabarty and Albert 2011). Central America is composed of three principal geological terrains; Nuclear Central America in the region of modern Honduras, El Salvador and Nicaragua dating from the Upper Cretaceous (100–66 Ma); Southern Central America in the region of modern Costa Rica dating from the Paleogene (66–23 Ma); and the Panamanian Isthmus in the region of modern Panama dating from the Late Miocene to Pliocene (12–3 Ma; Pindell and Kennan 2009; Bacon et al. 2015; Montes et al. 2015).

For the L + study we designed a dynamic model of Caribbean plate evolution using paleogeographic data summarized in Pindell and Kennan (2009). estimates were generated between five geographic areas (A–D) over three time intervals (Fig. 1, I–III). The five areas were delineated using published biogeographic criteria (Chakrabarty and Albert 2011; Matamoros et al. 2012, 2015): A) cis-Andean South America (Orinoco and Amazon basins); B) trans-Andean South America (Magdalena, Atrato, San Juan, and Baudo basins); C) Southern Central America (including Costa Rica and Panama); D) Nuclear Central America (Chortis Block including most of Nicaragua, Honduras, Guatemala, and El Salvador); and E) Greater Antilles (Cuba, Hispaniola, Puerto Rico, Jamaica). Species geographic distributions were estimated from museum collection records.

The three time intervals encompass geological events previously hypothesized to have affected vicariance and dispersal across the Caribbean plate. Time interval I encompasses the Upper Cretaceous to Eocene (80–35 Ma) movement of the Proto-Antilles between South America and Central America. Time interval II encompasses the Greater Antilles-Aves Ridge (GAARlandia) land bridge purported to have served as a transient terrestrial connection between northern South America and Greater Antilles at about the Eocene–Oligocene boundary (ca. 34–33 Ma; Iturralde-Vinent and MacPhee 1999). Time interval III encompasses the Middle Miocene to Pleistocene (12–3 Ma) rise of the Isthmus of Panama. In the paleogeographic model dispersal was disallowed between South America and the Western margin in Time intervals I and II (parameters .

Observed–Dispersal Rates

We define “dispersal rate” () as the number of biogeographic range expansion events per million years in taxa leaving descendants that persist to the Recent. We measure in units of “wallaces” (wa) to honor the contributions of Alfred Russel Wallace to the field of biogeography. The wa is introduced here as a standard metric to compare observed and expected rates of species-range expansion. This wa can also be used to compare quantitative dispersal rate estimates generated from biogeographic studies of different taxa across time and space, as for example, in the emerging fields of parametric biogeography (e.g., Landis et al. 2013; Peter Linder et al. 2013; Meseguer et al. 2014) and invasion biology (e.g., Schurr et al. 2012; Phillips 2012; Travis and Dytham 2012).

Here dispersal refers to the colonization of new areas outside an established species range, in contrast to dispersion, which refers to movements of organisms within an established species range (Armstrong 1977; Platnick 1979; Lomolino et al. 2010). Observed dispersal rate estimates ( were generated in L+ from 100 replicate maximum-likelihood searches, each initialized with a randomly produced set of dispersal rates (Table 1). Maximum-likelihood searches replicates were run in parallel on a 680 node shared Linux cluster operated by the Louisiana Optical Network Initiative until each search had converged. Results were pooled across replicates to find the parameter values with the best AIC value. Sensitivity analyses to investigate the effect of dispersal constraints on estimates of ancestral areas were also conducted in L+, employing dispersal rates selected from across a broad range of plausible values (102–10−10wa). Because L does not perform reliably with dispersal input values of 0, low-end values were input as 10−10wa.

Table 1.

Observed dispersal rates ( for freshwater fish clades across the Caribbean plate margins

Notes: Absolute values (in wa) for eight rate parameters (Greek letters as in Fig. 1) estimated by L+, reported as mean ± 1 SEM, median, and best ML values from 100 replicate runs. Global (tree-wide) extinction rates () and negative log likelihood (Log) values also reported for each clade.

Expected Dispersal Rates

Expected dispersal rates ( were generated from nine alternative paleogeographic models, employing two parameters based on island biogeography theory (dispersal rates proportional to target area, and inversely proportional to distance), and one parameter from freshwater fish ecology (dispersal rates proportional to water discharge volume of riverine source) (Matthews 1998). These paleogeographic models parameterize two properties (water discharge volume and river mouth location) of the Proto–Orinoco–Amazon (POA) paleoriver system that drained much of northern South America into the Caribbean Sea throughout most of the Cenozoic, before the Late Miocene (ca. 10 Ma) formation of the modern Orinoco and Amazon basins (Wesselingh and Hoorn 2011). We generated values from the equation: where is the relative land areas of the Eastern or Western plate margins, taken as proportions of modern values, is the relative water discharge volume of the POA paleobasin, also taken as a proportion of modern values, is the average paleogeographic distance from the POA river mouth to the closest point of the Eastern or Western plate margins during the time interval, and is a scaling constant in units of events/m3.

Values of were held constant in all analyses under the assumption that dispersal probability is a linear function of dispersal-target area and water discharge volume. This assumption can be relaxed in future implementations. Relative values of Eastern and Western margin areas were taken as the amount of emergent land-area of the Caribbean plate margins (Iturralde-Vinent and MacPhee 1999), under two different model conditions. The first condition is constant area (, assuming a single area value for all three time intervals. The second condition is variable areas ( assuming different areas for the three time intervals with a two-fold change in area between each time interval, associated with hypothesized sea-level shifts at the Eocene–Oligocene boundary (34 Ma) and from the Lower to Middle Miocene (12 Ma). Values of the POA river discharge volume were also taken as a percent of the combined modern Amazon + Orinoco discharge, assuming constant ( or variable ( river discharges over the three time intervals assessed. was modeled as two-fold changes in POA discharge associated with major orogenic events in the northern Andes through the Cenozoic (deGamero 1996; Gregory-Wodzicki 2000; Mora et al. 2010; Albert et al. 2011; Escalona and Mann 2011). The orogenic events are: Paleocene–Eocene at 59–55 Ma, Incaic at 42–35 Ma, Quechua 1 at 23–17 Ma, Quechua 2 at 9.5–7.0 Ma, and Quechua 3+4 at 6–0 Ma (Albert and Reis 2011, fig. 1.7). Paleogeographic distances were estimated directly from paleogeographic reconstructions in Pindell and Kennan (2009).

Results

Phylogenetic trees and ancestral area estimates are summarized in Figure 2 and Supplementary Figure S1, available on Dryad at http://dx.doi.org/10.5061/dryad.11cv0. Results for cichlids depict interrelationships for the 97 heroine species of the total of 172 cichlid species used in the analysis. The results of the L+ ancestral area estimates reveal a geologically persistent dispersal corridor along the Eastern Margin of the Caribbean plate during the Paleogene (66–23 Ma). The earliest time frame for dispersal to the Antilles occurred from 48.9 to 42.6 Ma in Heroini, from 59.6 to 49.5 Ma in Poeciliinae (Fig. 3). Ancestral areas estimated for heroine and poeciliine fishes using 12 alternative models of biogeographic range evolution in B (Matzke 2013a, 2013b) are reported in Supplementary Figure S2, available on Dryad at http://dx.doi.org/10.5061/dryad.11cv0, and likelihood results and AIC statistics from these analyses in Table 2. Quantitatively, the best-fit analyses for Heroini employed the DECj model of biogeographic range evolution, and the B-A-L set of parameter treatments, and the best-fit analysis for Poeciliinae employed the DECj model of biogeographic range evolution, and the D-L set of parameter treatments. Qualitatively, the L+ ancestral area estimates for tree branches in both Heroini and Poeciliinae (Fig. 3) at transitions among the five geographic areas (Fig. 1) are identical to those recovered by all 12 alternative models of biogeographic range evolution (Supplementary Fig. S2, available on Dryad at http://dx.doi.org/10.5061/dryad.11cv0).

Figure 2.

Normalized observed dispersal rates ( of Heroini (diamond symbols and blue lines) and Poeciliinae (square symbols and red lines) through time and across the Caribbean plate margins. estimates from L+ normalized to maximum observed values of each rate parameter (Greek letters as in Fig. 1). Note estimates along the Eastern margin (, indicate high connectivity between northern South America and Nuclear Central America from the Cenomanian to Eocene (ca. 80–35 Ma), but are not consistent with GAARlandia hypothesis of high dispersal rates at or around the Eocene–Oligocene boundary (ca. 34–33 Ma). Note also that estimates for the Western margin (, indicate high connectivity only during the Neogene (12–0 Ma), associated with biotic interchanges across the emerging Isthmus of Panama. These results illustrate a pattern of coordinate dispersal, that is, qualitatively similar values between the two fish clades in 11 of the 12 rate parameters across the Caribbean plate margins, differing substantially only at .

Figure 3.

L+ ancestral area estimates reveal a geologically persistent dispersal corridor along the Eastern margin of the Caribbean plate. Time-calibrated phylogenies of heroine cichlids (left) and poeciliine live-bearers (right) from beast analyses newly reported herein. Vertical gray bars indicate earliest time interval for the colonization of Nuclear Central America via the Eastern plate margin. Colors of branches and ancestral ranges indicate geographic areas (A–E) of Figure 1; biotic dispersals indicated by color changes along branches. Squares at branch corners indicate most likely ancestral ranges, of either one (i.e., A) or multiple (i.e., ABDE) areas.

Table 2.

Summary of likelihood results and AIC statistics for 12 alternative models of area-evolution

Note: (A) Heroini and (B) Poeciliinae. All analyses performed using the program B (Matzke 2013a, 2013b).

estimates for Eastern margin parameters ( indicate high connectivity between northern South America and Nuclear Central America during Time interval I (80–35 Ma), but low connectivity during Time interval II (34–13 Ma) and Time interval III (12–0 Ma). The estimates for Western margin parameters (, indicate high connectivity among trans-Andean, Southern Central American, and Nuclear Central American areas during the Neogene (12–0 Ma; Time Interval III), suggesting biotic interchanges during the Plio-Pleistocene across the Isthmus of Panama (Bermingham et al. 1997; Reeves and Bermingham 2006; McCafferty et al. 2012; Alda et al. 2013). On the other hand, estimates of for the Western margin indicate more dispersal among trans-Andean (area B), Panama (area C), and Nuclear Central America (area D) during the Neogene (Time III: 12–0 Ma). The relatively high values of and in Figure 3 indicate successful colonization events from Nuclear Central America (area D) and Southern Central America (area C) to trans-Andean northern South America (area B), events referred to as “Isthmian biogeographic reversals” by Chakrabarty and Albert (2011).

Correlations between estimates of and were assessed as coefficients of determination ( values from linear regressions) for each of the three time intervals of Figure 1, with values multiplied by , the duration of the time interval (in MY). Adjusted ( bar-squared) values were used to adjust for different numbers of model parameters (degrees of freedom) in correlations of data pooled for Heroini and Poeciliinae. The nine alternative paleogeographic models (Table 3) predicted relatively high expected dispersal rates ( along the Eastern margin during the late Upper Cretaceous and early Paleogene, and relatively low values along the Western margin during this same time interval. The variance of estimates increases with age due to accumulating uncertainties in relative POA discharges.

Table 3.

Expected dispersal rates ( for freshwater taxa across the Caribbean plate margins

Notes: (A) Absolute values (in wa) from nine alternative paleogeographic models. values from equation: (/, where is the relative land area of Eastern or Western margins, is the relative water discharge volume of Paleo–Orinoco–Amazon river, represents the proportion of modern values, is the closest paleogeographic distance between the coastlines of northern South America and plate margins, and is a scaling constant (see Methods for details). (B) values normalized to 1.0 from the maximum value ( of each model across a given plate margin. (C) Paleogeographically derived data used in model parameterization.

Geographic distance () is the single best predictor of among the paleogeographic parameters examined in this analysis, with values approximating unity in regressions against all four dispersal rate parameters (– in both Heroini and Poeciliinae (Table 4). Two other paleogeographic dispersal parameters, (but not and (but not , significantly predict dispersal rates on the Eastern margin (, . These results imply that the land surface area of the Greater Antilles and the relative volume of the POA discharge strongly affected the probability of overseas colonization. The lack of an association between and area variability ( on the Eastern margin may result from the approximately perpendicular orientation of the plate margin arc to the north coast of South America. In this configuration, even large changes in total emergent land area do not substantively alter the target area for overseas dispersal. However, several paleo land arcs potentially shortened this distance between the American continents and the Antilles (or even closed it) at several times since the Cretaceous (33). The correlations between and we recovered for Heroini are almost identical to correlations in these dispersal rates we obtained from a separate L+ analysis of a previously published phylogenetic analysis of Heroini (Říčan et al. 2013). The main differences between the two analyses are in the significance values for pooled data in model parameters and (Table 5).

Table 4.

Coefficients of determination () of regressions between normalized values of and for freshwater fish clades across the Caribbean plate margins

Notes: values from Table 1 and values from Table 3. Dispersal rate parameters (Greek letters) as in Figure 1. values for Heroini and Poeciliinae >0.950 in bold type; adjusted values for pooled data in bold type are significant at ( bar-squared >0.533, = 6). Paleogeographic distance () is the best predictor of on both plate margins, and that time-constant areal extent (in km2) of dispersal target area (, and time-varying volume of the POA river discharge (, are good predictors of on the Eastern margin, but not the Western margin. These results are consistent with overseas dispersal via a freshwater plume to the Antilles during the Paleogene, but not to Central America during the Neogene. These regression results are almost identical for Heroini and Poeciliinae, suggesting similar paleogeographic factors resulted in coordinated dispersal from South America to the Greater Antilles and Central America.

Table 5.

Coefficients of determination () between and for Heroini estimated from the phylogenetic hypothesis of Říčan et al. (2013)

Notes: Correlations differ from those in Table 4 only in the significance values for pooled data in models parameters and , with slightly lower correlation coefficients for , and slightly higher correlations for .

On the Western margin, paleogeographic distance () is the only significant predictor of dispersal rate (parameters , . Neither or alone or in combination are good predictors. Models combining distance, area, and discharge volume do not provide substantially better fits than the single-parameter models. The qualitative fits of alternative paleogeographic dispersal models to parameters are identical in the L+ analyses of Heroini and Poeciliinae; that is, these clades exhibit similar patterns of paleogeographic connectivity through time.

Discussion

Our main finding is that paleogeographic distance ( is the strongest predictor of observed dispersal rates ( in heroine cichlids and poeciliine live-bearers (Table 4). This result is consistent with overseas dispersal by means of rafting on a freshwater plume to the Antilles during the Paleogene, but not to Central America during the Neogene. Due to the westward movement of South America relative to the Caribbean plate, the mouth of the POA river reached its closest proximity to the Eastern (Antillean) margin during the Paleogene. Schubert (1935) and other early investigators (Forster 1778; Matthew 1919; Darlington 1938) recognized the importance of geographic distance limiting the dispersal of continental organisms to oceanic islands, and distance from source species pools is a principal parameter of the theory of Island Biogeography (MacArthur and Wilson 1967; Simberloff 1974; Whittaker and Fernández-Palacios 2007). The expectation that dispersal probability declines with the inverse square of geographic distance is predicted by ecological theory (Harper 1977; Viswanathan et al. 1999; Shaw et al. 2006; Cushman and Landguth 2010). An inverse square relationship of ecological dispersion with distance is widely reported in plants and diseases (Willson 1993; Turnbull et al. 2000; Mundt et al. 2009), and to a lesser extent in animals (Griffin 1952; Murray and Gillibrand 2006; Meats and Edgerton 2008). The inverse square law has only occasionally been applied in the context of island biogeography (Quinn et al. 1987; Gardner and Engelhardt 2008).

[T]he chances of dispersal of organisms across a barrier vary inversely as the square of distance (Schubert 1935).

The probability of rare events increases with longer time frames (Landres et al. 1999).

Our results also show a strong positive relationship between overseas dispersal probability in heroines and poeciliines, and the volume of the POA freshwater plume (Table 4). This finding supports the hypothesis that paleo-river discharges influence the formation of dispersal-assembled regional freshwater fish faunas in the Neotropics (Albert et al. 2006). The modern Amazon river has the largest freshwater discharge in the world, with an annual mean discharge rate of about 193,00 m3 s−1, an annual volume of , and a maximum areal extent of (Milliman and Meade 1983; Grodsky et al. 2014). The modern Orinoco has the fourth largest annual discharge at 31,00 m3 s−1 (Perry et al. 1996). The seasonally recurrent plumes of these large Neotropical rivers form large offshore lenses of low-salinity water along the northern margin of South America, varying from 80 to 200 km wide and 3 to 10 m thick (Muller-Karger et al. 1988). Direct measurements of sea surface salinities from buoys and satellite remote sensing show sea water with riverine values (absorption coefficients close to 1.0) about 580 km from the nearest point to the South American coastline (Hu et al. 2004). By comparison the south shore of Hispaniola is about 585 km from the nearest points on the South American coastline. The Amazon freshwater plume also has strong effects on the ecology (Johns et al. 2014) and biogeography (Ludt and Rocha 2015) of Caribbean marine fish species. Overseas fish dispersal events may also have been facilitated by ephemeral island arcs that shortened the distance between landmasses, and the prevailing southeast-to-northwest direction of the circumtropical paleocurrent and Western Atlantic hurricane track.

The results of the model-fitting (Table 4) imply a history of overseas dispersal events across the Caribbean plate margins, coordinated and constrained by the paleogeographic configuration of landmasses, and by prevailing and perennial vectors of atmospheric and oceanic circulation operating over geological time intervals. Coordinated dispersal refers to the influence of similar paleogeographic constraints on the distribution of taxa in space, and does not necessarily imply coincident dispersal events in time (Lieberman and Eldredge 1996; Albert and Carvalho 2011). In this regard, coordinated dispersal differs from range vicariance as alternative explanations for concordant biogeographic patterns, because vicariance but not coordinated dispersal may be falsified by evidence for the pseudocongruent timing of events (Cunningham and Collins 1994; Taylor et al. 1998; Donoghue and Moore 2003).

The timing of dispersal events recovered here supports the hypothesis that cichlids and poeciliids became established in Central America during the early Cenozoic. Because most GACA freshwater fish clades have closest relatives in South America, early workers speculated these taxa originated from dispersal over the Panamanian land bridge about 3.0–3.5 Ma (Lomolino et al. 2010). Yet it has long been hypothesized that some obligate freshwater fish taxa, cichlids and poeciliids in particular, may have dispersed to Central America before the rise of the Panamanian Isthmus (Myers 1949, 1966; Miller 1966; Bussing 1976, 1985) by overseas dispersal. Recent phylogenetic studies of Neotropical freshwater fishes have supported a pre-Isthmian date for expansion into Central America (Smith and Bermingham 2005; Chakrabarty 2006; Hrbek et al. 2006; Chakrabarty and Albert 2011; Říčan et al. 2013; Matamoros et al. 2015).

In terms of spatial relationships, our results are consistent with the general area relationships, but not the precise timing, predicted by the GAARlandia hypothesis (Iturralde-Vinent and MacPhee 1999). Our results place the earliest colonizations of Central America in the Paleocene–Eocene (ca. 58–45 Ma), long before the posited age of the GAARlandia land bridge (ca. 34–33 Ma). These results suggest that GAARlandia, if it ever existed as a subaerial land bridge (Ali 2012), was a transient structure embedded within the larger geological history of the Eastern margin. These results highlight the role of the Eastern margin as a dispersal corridor for freshwater fishes between portions northern South America, the Greater Antilles, and the Maya (Yucatán) block, and also possibly subaerial portions of the Aves Ridge as island arcs at various times in the early Cenozoic and of varying lengths (Iturralde-Vinent and MacPhee 1999; Iturralde-Vinent 2006). These results are consistent with recent discoveries demonstrating the role of improbable overseas journeys by means of rafting on the formation of insular biotas (Measey et al. 2007; De Queiroz 2014).

Methods for estimating macroevolutionary rate parameters from molecular phylogenies are still in an active stage of development and the precise biological meaning of the rate estimates reported here is incompletely understood (Rabosky 2010; Morlon et al. 2011; Stadler 2011; Pyron and Burbrink 2013). For example, the extinction rates reported here for freshwater fishes in the northern Neotropics (∼1E–2; Table 1)are lower than extinction rates reported for other fish taxa based on the fossil record (McShea 1994; Ricciardi and Rasmussen 1999). As in all parametric methods in macroevolution, L+ estimates of and are based on available phylogenetic and biogeographic data, and not the complete sets of (living and extinct) taxa and geographic distributions that characterize the full evolutionary history of a group (Chakrabarty 2004; Hrbek et al. 2006). These rate estimates are minimally sufficient to explain current phylogenetic and biogeographic distributions of tip species. Despite these important filters, parametric biogeographic analysis can reveal important patterns in organismal diversification in time and space (Ree and Sanmartin 2009; Pyron and Burbrink 2013; Pennell et al. 2014).

L+ dispersal rate estimates reported here (Table 1) were calculated in the context of a formal paleogeographic landscape model constructed using information on the timing and location of geological events (Pindell and Kennan 2009). This geological knowledge can strongly affect estimates of some biogeographic rate processes. For example, knowledge that the Isthmus of Panama arose during the late Neogene (12–3 Ma) constrains numerical estimates of dispersal rates between Nuclear Central America and South America during Time interval II, such that values of from the unconstrained analysis are relatively high (3.19e10–2) when compared with the value of 0.00 in the paleogeographically constrained analysis. L+ estimates from paleogeographically constrained and unconstrained analyses were also qualitatively different for three rate parameters (, on the Eastern Margin during Time intervals II and III (Table 6).

Table 6.

L+ estimates of for Heroini and Poeciliinae across the Caribbean plate margin, with and without paleogeographic constraints

	Paleogeographic model
Parameter	Constrained	Unconstrained	Constrained	Unconstrained
	Heroini		Poeciliinae
Eastern margin
α_1	1.37E−01	1.00E+02	3.49E−01	1.28E−01
α_2	3.15E−03	3.58E−28	2.66E−03	3.02E−40
α_3	1.04E−32	2.43E−25	4.33E−03	5.28E−03
β_1	1.19E−02	6.76E−03	1.51E−01	1.25E−02
β_2	2.72E−03	2.44E−30	7.57E−31	3.74E−23
β_3	3.92E−03	5.57E−29	9.31E−36	2.75E−21
Western margin
δ_1	0.00E+00	3.83E−16	0.00E+00	1.69E−26
δ_2	0.00E+00	1.45E−28	0.00E+00	3.82E−23
δ_3	1.07E−02	3.27E−02	5.52E−02	4.85E−02
ε_1	0.00E+00	7.71E−27	0.00E+00	4.73E−25
ε_2	0.00E+00	3.19E−02	0.00E+00	1.30E−02
ε_3	4.94E−03	1.92E−03	8.59E−03	2.83E−03
Ensemble
μ	1.32E−02	1.39E−02	1.92E−02	1.25E−02
nLog L	9.58E+01	1.05E+02	1.32E+02	1.26E+02

Notes: Paleogeographically constrained values of Western margin rate parameters at time intervals I and II (, , , ) set to 0.00. Unconstrained rate values estimated from the time-calibrated trees and extant species distributions alone (Fig. 1). ML values from 100 replicate runs. Ensemble (tree-wide) extinction rates () and negative log likelihood (Log ) values also reported. Rates estimates from the constrained and unconstrained analyses are significantly correlated (Heroini, , ; Poeciliinae, , ). However, rates estimates for the constrained and unconstrained analyses differ by many orders of magnitude (values in bold) for three rate parameters on the Eastern margin (, ), and three rate parameters on the Western margin (, , ).

Phylogenetic and biogeographic data reviewed here suggest the Pleistocene Isthmian exchange was of limited taxonomic scope for freshwater fishes, bidirectional, with taxa moving in both directions, and asymmetrical, with a preponderance of taxa moving south into the Pacific slope of Colombia (Chakrabarty and Albert 2011). The rise of the Isthmus is best viewed as the most recent of many geological and geographic events involved in the formation of the modern Central American and trans-Andean ichthyofaunas. The Isthmus is only one piece of a richly complex puzzle that is the biogeographic history of this region.

Author's Contributions

V.A.T., S.M.D.-S., and J.S.A., conceived the ideas; V.A.T., P.C., and W.A.M. collected the data; S.M.D.-S. wrote the computer code; all authors analyzed the data; J.S.A. led the writing.

Supplementary Material

Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.11cv0.

Funding

This work was supported by the United States National Science Foundation [DEB 0916695 and 1354511 to P.C., and DEB 0614334, 0741450, and 1354511 to J.S.A.]; the Fundação de Amparo à Pesquisa do Estado de São Paulo [2012/09990-0 to V.A.T.]; and the Louisiana Education Quality Support Fund [2011-14-RD-A-27 to S.M.D.-S.].