- Research article
- Open Access
Paleoclimatic modeling and phylogeography of least killifish, Heterandria formosa: insights into Pleistocene expansion-contraction dynamics and evolutionary history of North American Coastal Plain freshwater biota
© Bagley et al.; licensee BioMed Central Ltd. 2013
- Received: 2 September 2012
- Accepted: 13 September 2013
- Published: 9 October 2013
Climatic and sea-level fluctuations throughout the last Pleistocene glacial cycle (~130-0 ka) profoundly influenced present-day distributions and genetic diversity of Northern Hemisphere biotas by forcing range contractions in many species during the glacial advance and allowing expansion following glacial retreat ('expansion-contraction’ model). Evidence for such range dynamics and refugia in the unglaciated Gulf-Atlantic Coastal Plain stems largely from terrestrial species, and aquatic species Pleistocene responses remain relatively uninvestigated. Heterandria formosa, a wide-ranging regional endemic, presents an ideal system to test the expansion-contraction model within this biota. By integrating ecological niche modeling and phylogeography, we infer the Pleistocene history of this livebearing fish (Poeciliidae) and test for several predicted distributional and genetic effects of the last glaciation.
Paleoclimatic models predicted range contraction to a single southwest Florida peninsula refugium during the Last Glacial Maximum, followed by northward expansion. We inferred spatial-population subdivision into four groups that reflect genetic barriers outside this refuge. Several other features of the genetic data were consistent with predictions derived from an expansion-contraction model: limited intraspecific divergence (e.g. mean mtDNA p-distance = 0.66%); a pattern of mtDNA diversity (mean Hd = 0.934; mean π = 0.007) consistent with rapid, recent population expansion; a lack of mtDNA isolation-by-distance; and clinal variation in allozyme diversity with higher diversity at lower latitudes near the predicted refugium. Statistical tests of mismatch distributions and coalescent simulations of the gene tree lent greater support to a scenario of post-glacial expansion and diversification from a single refugium than to any other model examined (e.g. multiple-refugia scenarios).
Congruent results from diverse data indicate H. formosa fits the classic Pleistocene expansion-contraction model, even as the genetic data suggest additional ecological influences on population structure. While evidence for Plio-Pleistocene Gulf Coast vicariance is well described for many freshwater species presently codistributed with H. formosa, this species demography and diversification departs notably from this pattern. Species-specific expansion-contraction dynamics may therefore have figured more prominently in shaping Coastal Plain evolutionary history than previously thought. Our findings bolster growing appreciation for the complexity of phylogeographical structuring within North America’s southern refugia, including responses of Coastal Plain freshwater biota to Pleistocene climatic fluctuations.
- Coastal Plain
- Mismatch Distribution
- Incomplete Lineage Sorting
- Coalescent Simulation
- Atlantic Coastal Plain
Present-day distributions of many of Earth’s biotas reflect the profound influence of climatic and sea-level fluctuations during the glacial-interglacial cycles of the Pleistocene (2.58-0.01 million years ago, Ma; ) [2–4]. Notably, eight glacial advances since 740 ka have produced more extreme 100,000-year environmental changes than those preceding 1 Ma . Because glacial-interglacial transitions occurred in an evolutionary blink of an eye, many extant species have spent most of their recent evolutionary histories under glacial rather than short-lived (~10,000-20,000 years) interglacial conditions . Therefore, extreme environmental changes associated with the Last Glacial Maximum (LGM; ~22-19 ka) may have been the most recent climatic factors that affected present-day patterns of biodiversity prior to the global proliferation of humans and anthropogenic environmental disturbance [2, 6].
The biogeographical consequences of Pleistocene glaciations are generally thought of in terms of direct effects of glacial cover and temperature change on Northern Hemisphere biotas. This is because over 80% of glacial ice on Earth resided in the Northern Hemisphere during the LGM, resulting in the destruction and reorganization of continental landscapes in both the Old and New World e.g. [2, 7]. While few species from glaciated areas persisted in situ by adapting to changing local environments, many were actively displaced (tracking suitable habitat), passively displaced (local extinction, regional extirpation), or went extinct [3, 8]. Indeed, abundant paleontological and palynological evidence indicates that many terrestrial species of North America and Europe contracted their ranges during the LGM to one or more lower-latitude (or altitude) refugia, then expanded northward to recolonize much of their present-day ranges during the warmer Holocene interglacial [9–14]. This scenario is known as the 'expansion-contraction’ model of Pleistocene biogeography and has reached paradigm status in biology .
Many patterns of genetic variation are consistent with predictions from a model of range contraction during glacial advance and rapid post-glacial expansion. For example, phylogeographical analyses across European flora and fauna have identified the Iberian, Italian and Balkan peninsulas as three main southern refugia from which postglacial recolonization of northern Europe took place [9, 16, 17]. In more recent work, the tools of phylogeography and ecological niche modeling e.g.  have been joined to refine models of Pleistocene expansion and contraction. Such analyses have revealed complicated patterns of northern and southern LGM refugia and complex histories of speciation, dispersal, secondary contact, and demographic expansion of species within the southern refugia themselves, particularly in Europe [15, 18–24].
The Gulf-Atlantic Coastal Plain of eastern North America presents an ideal system in which to study interactions between geology, climate and diversification of a unique biota. This unglaciated region has figured prominently in phylogeographic research, and has revealed several major genetic breaks and confirmed local glacial refugia for northern taxa and endemic coastal species [29, 30]. Many Coastal Plain species possess evolutionary lineages that have diverged in an east-to-west direction across common biogeographic barriers, including peninsular Florida which separates the maritime Atlantic-vs.-Gulf coast biota, the Apalachicola River (Gulf-Atlantic slope drainages), Mobile Bay, and the Mississippi River [30–32]. Genetically identified breaks in this region also coincide with species distributional breaks  and thus likely reflect historical vicariance events, or recurrent vicariance coupled with extinction-recolonization of populations [29, 31, 34]. However, phylogeographical studies also indicate that species persisted through the LGM in refugia on either side of the Apalachicola basin in Florida and Texas-western Louisiana [29, 31]. Between these regions, there are north-south-trending zones of hybridization and secondary contact, including one of Remington’s 'suture zones’ between Alabama and northwest Florida, indicating a meeting place for lineages undergoing northward postglacial expansion [30, 35]. As is the case for the patterns in evidence for expansion-contraction dynamics in unglaciated regions of Europe [2, 4, 9], evidence for these refugia in North America has been developed mainly from terrestrial taxa, e.g. mints (e.g. Conradina[36, 37]) and yellow-poplar (Liriodendron tulipifera). These terrestrial data suggest that other codistributed species may have experienced similar range dynamics fitting the expansion-contraction model, possibly involving similar refugia. However, while regional genetic breaks are well documented, Pleistocene evolutionary responses of Coastal Plain aquatic biota remain relatively uninvestigated.
Here, we focus on least killifish, Heterandria formosa Girard 1859, a species which presents several advantages as a historical biogeographic study system that likely exhibited distinct response(s) to the climatic upheavals of the Late Pleistocene. First, these small (12-30 mm) livebearing freshwater fish (family Poeciliidae) are restricted to low-elevation Coastal Plain areas of subtropical humid climate zones (hot, wet summers and mild winters), and the tropical Everglades . This suggests climate is likely an important factor limiting H. formosa distribution. This species is also intolerant to cold-seasonal temperatures, making it an ideal candidate for ecological niche modeling and for testing the expansion-contraction model. Second, H. formosa span the Florida peninsula, an important glacial-stage refugium [29, 35]; thus, it is plausible that populations were displaced to warmer south-peninsula areas during the LGM. Third, H. formosa display ecological characteristics consistent with the potential for rapid population expansion over evolutionary timescales, including (i) rapid reproduction and short generation time (T ≈ 0.33 yr ); (ii) a range of tolerances to different abiotic conditions [39, 40], and (iii) female capacity for livebearing and sperm storage, suggesting individuals or groups of migrating females could found populations during rapid post-LGM expansions . Fourth, H. formosa is endemic to the Gulf-Atlantic Coastal Plain and thus also offer the opportunity to test the generality of the standard regional vicariance model and phylogeographical breaks documented in other taxa.
A previous rangewide study by Baer  using allozymes inferred genetic barriers in H. formosa between the Western Coastal Plain (WCP) and Atlantic Coastal Plain (ACP) regions (Figure 1), around the Suwannee River, but not between north Florida and the ACP. Another break occurred between Gulf-draining Waccasassa and Withlacoochee Rivers, creating a clade of Louisiana plus south Florida samples. Baer  also hypothesized an important role for population expansion and cross-peninsula gene flow in influencing current levels of genetic variation—including recent founding of ACP populations after Early Pleistocene high seas (≤1.5 Ma; 'northeast colonization’ hypothesis), supported by lower genetic diversity and a lack of isolation-by-distance only among ACP populations. Yet while Baer’s  study presents a broad outline for understanding H. formosa historical biogeography, the limited levels of genetic variation he found and the more modest methods available at the time limited the hypotheses that could be tested. The current analytical framework in historical biogeography, including our ability to combine ecological niche modeling (e.g. to infer Pleistocene-Holocene distributions and generate spatially explicit hypotheses) with coalescent simulations (e.g. used to discriminate among historical scenarios [43, 44]), provides a more rigorous framework for biogeographical inference in this species [45–47].
In this study, we integrate ecological niche modeling with paleoclimatic data, phylogeography, and coalescent simulations to infer the historical biogeography of H. formosa. Specifically, we use H. formosa to test the hypothesis that members of the Coastal Plain freshwater biota fit predictions of the expansion-contraction model (as opposed to southern crossroads, or regional vicariance and northeast colonization hypotheses) by testing for predicted patterns of Late Pleistocene range shifts, refugia and recolonization patterns and their spatial-genetic effects. In addition to analyzing new mtDNA and nDNA sequences, we re-analyze available allozyme data , which provide a novel temporal spectrum and improve our tests of genetic diversity predictions, including latitudinal patterns not previously investigated. We show that these new data, combined with paleoclimatic modeling in an integrative approach, support a hypothesis of Pleistocene range expansion-contraction dynamics in least killifish.
Sampling and laboratory methods
Heterandria formosa collection sites, sample sizes and genetic diversity measures
Population (specimen no. code†)
Latitude (° N)
Longitude (° W)
Bayou Nezpique (BNE x LA)
Bayou Nezpique, south of I-10, ~30 km east of Lake Charles, Cameron Par., LA
Lake Pontchartrain (PON x LA)
Lake Pontchartrain, LA
Bogue Chitto (BOG x LA)
Bogue Chitto R. overflow pool along LA Hwy 21 south of Sun, St. Tammany Par., LA
Porter's River (POR x LA)
Porter's R. at Davis Landing Rd., St. Tammany Par., LA
Wolf Creek Swamp (WCMS x )
Wolf Cr. swamp at Belle Ferry Rd., Harrison Co., MS
Alabama River trib. (ALA x AL)
Unnamed trib. to Alabama R. (Bailey Cr. origin), Wilcox Co., AL
Five Points (FIV x FL)
Unnamed trib. at Five Points, Tate's Hell State Forest, Franklin Co., FL
Crooked River (CRO x FL)
Crooked R. at Tate's Hell State Forest, Franklin Co., FL
Ochlockonee River (OCH x FL)
Ochlockonee R. at Blountstown Hwy, Leon Co., FL
Womack Creek (WOM x FL)
Womack Cr. at Apalachicola National Forest, Franklin Co., FL
1, 3, 11-14
Moore Lake (MOO x FL)
Moore Lake at Apalachicola National Forest, Leon Co., FL
Hill Swale (HIL x FL)
Hill Swale, just north of Ochlockonee Bay, Wakulla Co., FL
Trout Pond (TRO x FL)
Trout Pond at Apalachicola National Forest, Leon Co., FL
Cessna Pond (CES x FL)
Cessna Pond at Apalachicola National Forest, Leon Co., FL
Little Lake (LIT x FL)
Little Lake Jackson, west of Hwy 27, northwest Tallahassee, Leon Co., FL
Wakulla Springs (WAK x FL)
Wakulla Springs at Edward Ball Wakulla Springs State Park, Wakulla Co., FL
7, 9-10, 20
Lake Iamonia (IAM x FL)
Lake Iamonia at CR 155, Leon Co., FL
Shepherd Spring (SHE x FL)
Shepherd Spring at St. Marks National Wildlife Refuge, Wakulla Co., FL
McBride Slough (MCB x FL)
McBride Slough at Bloxham Cutoff Rd., Wakulla Co., FL
Lake Overstreet (LOV x FL)
Lake Overstreet at Alfred B. Maclay Gardens, Leon Co., FL
Newport Sulphur Spring (NEW x FL)
Newport Sulphur Spring, Wakulla Co., FL
2, 9, 15, 18
Natural Bridge (NAT x FL)
Natural Bridge at Natural Bridge Rd., Wakulla Co., FL
Gambo Bayou (GAM x FL)
Gambo Bayou at Lighthouse Rd. west of Stony Bayou Pool, Wakulla Co., FL
Tram Road (TRA x FL)
Tram Rd., 2 km east of Plum Orchard, Wakulla Co., FL
Wacissa River (WAC x FL)
Wacissa R., just north of Aucilla Wildlife Management Area, Jefferson Co., FL
Buggs Creek borrow pit (BUG x FL)
Buggs Creek barrow pit, FL
Wolf Creek (WOL x FL)
Wolf Creek at US 90, Jefferson Co., FL
Mckey Park (MCK x GA)
Mckey Park, Loundes Co., GA
Bevil Creek (BEV x GA)
Bevil Creek at Browns pond Loch Laurel Rd., GA
Robinson Creek (ROB x FL)
Robinson Cr. at CR 131, Columbia Co., FL
Ichetucknee blue hole (ICH x FL)
Blue hole, Ichetucknee State Park, Columbia Co., FL
Hillsborough River (HIR x FL)
Hillsborough R. at Tampa, Hillsborough Co., FL
River Styx (RIS x FL)
River Styx at CR 346, Alachua Co., FL
Saint Johns River (SJR2 x )
Unnamed trib. of the St. Johns R., FL
Newnan's Lake (NEWn x FL)
Newnan's Lake at southeast 16th Ave. boat ramp, Alachua Co., FL
Everglades (EVE x FL)
Everglades drainage, Dade Co., FL
Bahama Swamp (BAH x SC)
Bahama Swamp off Kato Bay Rd., Jasper Co., SC
Edisto River (CRO x SC)
Edisto River, SC
Back River (BAC x SC)
Back R. at SSR 503, east of Goose Creek, edge of US Navy station, SC
Cooper River trib. (WAD1 x SC)
Trib. to Cooper River at Wadboo, Berkeley Co., SC
Waccamaw River (WAR1 x SC)
Waccamaw River, SC
Trib. to Waccamaw River (WAR2 x SC)
Unnamed trib. to Waccamaw River, SC
Lumber River (LUM x SC)
Lumber River, SC
Cane Branch (CAN x SC)
Cane Branch on farm road at Francis Marion National Forest, Berkeley Co., SC
We isolated DNA using the Qiagen DNeasy96 tissue protocol (Qiagen Sciences, Maryland, USA). We amplified the mitochondrial cytochrome b (cytb) gene for each sample by PCR using Hrbek et al.’s  L14725 - H15982 primer pair. We also sequenced nuclear ribosomal protein S7 (RPS7) introns 1 and 2 and exon 2 for subsamples of 1-3 H. formosa individuals/site for up to 4 random sites per region within each clade (based on our multilocus gene tree; see Results), and Poeciliidae sp. samples. We amplified RPS7 using a nested PCR design, with primers 1 F - 3R.24 in the first reaction followed by internal primers 1 F.2 - 2R.67 and 2 F.2 - 3R in subsequent reactions . Primers 1 F and 3R are from Chow and Takeyama . Amplification conditions, sequencing reactions, clean-up, and sequence visualization followed Unmack et al. . We aligned sequences manually while viewing electropherograms in Sequencher™ 4.10 (Gene Codes Corp.). All sequences generated in this study were deposited in GenBank (accession numbers: KF632895-KF633114, cytb, and KF633115-KF633133, RPS7). Alignment lengths were 1140 nucleotide base pairs (bp) for cytb and 876 bp for RPS7. We supplemented these data with orthologous ingroup sequences from Lake Pontchartrain, Louisiana (N = 1) and Everglades, Florida (N = 1) retrieved from GenBank (Additional file 1: Table S1). Thus, our sampling encompassed H. formosa from 44 georeferenced sites (Table 1). In total, the final cytb database we analyzed consisted of 223 H. formosa and 3 Poeciliidae sp. sequences that we collapsed into unique haplotypes in DnaSP 5.10 . To create a multilocus alignment, we collated individuals with cytb and RPS7 sequences and joined them with analogous sequences from potential outgroups. Except multilocus and haplotype-based tests, analyses used the full cytb database (sometimes divided among groups); iteratively adding/removing samples with missing data (usually <5 bp) had no qualitative effect on results.
Ecological niche modeling
To test the prediction that H. formosa experienced a southward LGM range contraction and a northward range expansion from LGM to present, consistent with the expansion-contraction model, we conducted ecological niche modeling using the maximum entropy (Maxent) approach implemented in MAXENT 3.3.3 k . The Maxent method predicts locations of suitable habitat based on environmental-climatic characteristics and identifies abiotic and biotic factors that may limit species present-day distributions. We used Maxent because it exhibits high predictive performance .
We based our GIS-based niche-modeling analyses on 259 georeferenced H. formosa sites (Additional file 1: Figure S2) collated from our field collections (Table 1) and Global Biodiversity Information Facility (http://www.gbif.org) records. Sampling covered the species known distribution and was beyond sufficient given ≥10 sites permit accurate niche model construction . From WorldClim (http://www.worldclim.org; ), we assembled GIS coverages with a resolution of 30 arc-seconds (1 km2) for 19 bioclimatic environmental predictor variables (Additional file 1: Table S2) for present-day climates (1950-2000), plus analogous paleoclimatic data layers reconstructing past LGM (21 ka) and LIG (140-120 ka) environments. The LGM dataset (resampled from 2.5 arc-minutes resolution) was derived from the CCSM3 global circulation model of the Paleoclimate Modelling Intercomparison Project Phase II (PMIP2; ). The LIG dataset (originally 30-s resolution) was derived from climate simulations in Otto-Bliesner et al. . Layers were clipped to an extent of 24.5°N-35°N and 95°W-106°W prior to analyses. Unfortunately, high-resolution data on freshwater environments is not available for all three time-periods modeled, although such data would aid niche-based modeling of freshwater organisms. However, the CCSM3 data estimate environments over continental shelf areas exposed during the LGM; thus, we were able to reproject our full models over subaerial LGM landmass within the modeled area, which was critical for testing the expansion-contraction model. We ran niche models using different combinations of predictor variables to broadly explore the utility of different variables and their effects on model prediction. Ultimately, we included all 19 data layers in the final modeling procedures, because models showed no evidence of over-fitting and MAXENT is robust to correlations among variables (Additional file 1: Supplement S1).
We predicted the current (0 ka) geographical distribution of H. formosa, by generating an ecological niche model in MAXENT using the program’s default settings (e.g. 104 background points, logistic model with habitat suitability between 0 and 1), except we increased the number of iterations to 5000 to ensure convergence and we averaged results over 10 replicate crossvalidation runs. All data were used for model testing and training. We reprojected the resulting niche models on LGM and LIG layers and interpreted areas with high-predicted LGM bioclimatic suitability as most-likely refugial areas. Under the expansion-contraction model, we expected areas of predicted occurrence to be reduced during the LGM relative to LIG and present-day, restricted to a southern refugium. By contrast, refugia located on either side of the Apalachicola River would support vicariance-northeast colonization, and a pattern of range stability through time would support southern crossroads. Model evaluation was based on threshold-independent measures of performance, area under the Receiver Operating Characteristic curve (AUC) statistics; scores closer to 1 (maximum) indicate higher predictive ability and AUC > 0.5 indicates better-than-random model prediction . We generated response curves showing the impact of each variable alone on MAXENT prediction, and we studied 'multivariate environmental similarity surfaces’ (MESS) output by the program to assess extrapolation and effects of novel environments on prediction.
Population structure and genetic diversity
We evaluated DNA polymorphism levels by calculating segregating sites (S), h, Hd, and π in DnaSP for the full H. formosa cytb dataset, as well as populations and SAMOVA groups (see Results) that met a threshold-sampling criterion of N ≥ 8. Due to many samples not meeting our N ≥ 8 threshold-sampling criterion, mtDNA data were not used to test predicted genetic diversity patterns. Watterson’s estimator (θW; from S) of θ was calculated overall and for SAMOVA groups. Additionally, we calculated uncorrected pairwise mtDNA p-distances between individuals and groups in MEGA 5 .
We analyzed the mtDNA data to test for patterns of spatial-genetic subdivision predicted by the various phylogeographical hypotheses. The expansion-contraction model predicts genetic drift should have created distinct lineages subdivided between refugia—if multiple, allopatric refugia existed—especially over repeated glacial cycles [2, 15]. However, if northward recolonization occurred from a single southern refuge via 'leading-edge colonization,’ this should be reflected in a genetic distinction between the two areas, separating established (expanded) from blocked interior populations/lineages [2, 9]. If the recolonization wave passed through spatial constraints, e.g. coastlines, we expect to identify most recent common ancestors (MRCA) through the coalescent and genetic distinctions that mark the location of barriers encountered outside of a refuge . For what we term the 'vicariance-northeast colonization’ hypothesis (uniting Coastal Plain vicariance e.g. [29, 32] and northeast colonization  elements), we tested for genetic barriers in the Apalachicola River-region between ACP-FL and WCP samples, and between the Suwannee and Hillsborough Rivers along Florida’s Gulf Coast (our sampling analog to the area of Baer’s  Waccasassa-Withlacoochee break). However, we expected genetic barriers to be absent within the ACP . By contrast, a southern crossroads scenario would be supported by a pattern of no phylogeographical structure .
We tested these predictions using Monmonier’s  algorithm, implemented in BARRIER 2.2 , to identify genetic 'barriers’ across a Delaunay triangulation network overlaying the sampling sites, based on maximum Tamura–Nei genetic distances calculated in Arlequin 3.5  (1000 nonparametric permutations). To assess barrier support, we calculated bootstrap proportions from 100 bootstrapped barriers generated in BARRIER using 100 bootstrapped Tamura–Nei distance matrices made in PopTools ; bootstrap values >50 indicated strong support. Next, we used spatial analysis of molecular variance (SAMOVA), implemented in SAMOVA 1.0 , to define population clusters maximizing the proportion of total genetic variance due to differences among geographical groups (ΦCT). We performed SAMOVAs across K = 2-8 groups using pairwise differences, drawing from any of 100 initial random conditions. We independently tested the grouping schemes indicated by the SAMOVA and BARRIER analyses by conducting analyses of molecular variance (AMOVA) in Arlequin (1000 nonparametric permutations). Because BARRIER and SAMOVA are sensitive to small sample sizes, we pooled adjacent sites until an N ≥ 8 threshold-sampling criterion was reached per population, yielding a 22-population subset.
We also conducted additional spatial analyses. Using AMOVA, we tested for significant genetic partitioning among drainages. Here, we capitalized on our dense north Florida sampling and split populations across seven drainages, including the Ochlockonee (collections 7-15, 17, and 20), St. Marks (16, 19, and 21-22), Apalachee plus Goose Creek Bays (18, 23-24), Aucilla (25-27), Suwannee (28-31), Hillsborough (32), and St. Johns Rivers (34-35). We employed GENALEX 6.3  to examine predictions about patterns of isolation-by-distance, using Mantel tests  for correspondence between mtDNA genetic distances [FST/(1 - FST)] and geographic distances [ln (km)] among sites at rangewide and FL-region levels (N ≥ 8 threshold-sampling; 104 permutations; relationships confirmed using linear regression). We used straight-line geographic distances, not river distances, because several sites were bodies of water (e.g. lakes) that are isolated from any nearby rivers. Based on our niche models (see Results), we predicted much of the range would lack isolation-by-distance; however, under vicariance-northeast colonization, only ACP populations should lack isolation-by-distance .
We tested for genetic diversity patterning consistent with the expansion-contraction model using analyses of the allozyme dataset. Coalescent theory predicts that northward postglacial expansion from reduced Ne should have created genetic signatures of population expansion in recently founded (recolonized) populations. These include excesses of rare alleles and low-frequency mutations, reduced genetic diversity (homozygosity, shorter demographic histories) resulting in negative latitudinal clines, and a lack of isolation-by-distance (reviewed in [4, 9, 69, 70]). Conversely, long-term stable refugia should display isolation-by-distance and higher genetic diversity. We tested for predicted latitudinal clines in diversity by calculating expected heterozygosity (He; Nei’s expected heterozygosity) over all loci in GENALEX (from allele frequencies), then testing for a negative relationship with latitude [4, 9, 17] using linear regression models in PAST 2.02 . An absence of clinal diversity patterning would be more consistent with a southern crossroads scenario. We tested whether allozyme He and private allelic richness (hp) were greater within the niche modeling-inferred refugial populations (see Results) than throughout the putatively recolonized modern range using Mann-Whitney and Kolmogorov-Smirnov tests in PAST (due to non-normality, non-homogeneous variance). Using similar nonparametric tests, we tested whether WCP or ACP diversity was significantly lower than that of other populations combined, or populations within the putative refuge. And we tested for expected patterns of (lack of) isolation-by-distance using Mantel tests of unbiased Nei’s D genetic distances and ln geographic distances between populations, calculated in GENALEX at rangewide, WCP, ACP, FL, and within-refuge levels.
We tested for population expansions predicted by the expansion-contraction model (as per above) at levels of regional groups, SAMOVA groups (see Results) and local populations (N ≥ 8 threshold size) using complementary cytb analyses. We considered the expansion-contraction model 'strongly’ supported where at least two statistical tests supported expansion. First, we estimated Fu’s FS, and R2, and their 95% confidence intervals using coalescent simulations in DnaSP (testing significance with 104 replicates). To distinguish population expansion from purifying or positive selection, we tested each group, and the full cytb database, for mtDNA neutrality using McDonald and Kreitman’s  test and coalescent simulations of Fay and Wu’s H in DnaSP (104 replicates). We used Poeciliidae sp. as the McDonald–Kreitman test outgroup; other potential outgroups from Additional file 1: Table S1 yielded identical results (unpublished data). Second, we conducted mismatch distribution tests for regions/groups with no apparent subdivision (WCP, ACP and SAMOVA-group levels) in Arlequin (1000 parametric bootstrapping iterations) to see if pairwise nucleotide difference frequencies rejected sudden- or spatial-expansion models (given similar results, only sudden-expansion results are presented). We calculated Harpending’s raggedness index (r) as a test statistic measuring goodness-of-fit to the mismatch distributions. We calculated the time in generations since population expansion (t) using estimated expansion parameters and the equation τ = 2 μt (τ = coalescent time in mutations, μ = mutation rate/locus/year based on our Bayesian cytb rate estimate; see Results); t was converted to time in years using generation time. We expected mean t and/or 95% confidence intervals to overlap the LGM.
Phylogenetic relationships and coalescent-dating analyses
We expected multiple Pleistocene refugia to be supported by allopatric clades/unique networks [4, 15] with strong nodal/parsimony support, whereas unique clades or intra-network structuring between refugia inferred from niche modeling and recolonized areas would support leading-edge effects [2, 9]. Another phylogenetic prediction is that derived haplotypes (mutations at phylogeny/network tips) should become geographically localized outside of refugia with ancestral network populations found within refugia . The placement of one or more network roots outside of the inferred refugia would be consistent with a model of allele surfing  or extinction of trailing expansion edges . Here, southern crossroads is a null model with no obvious phylogeographical patterning within a single population-lineage.
We inferred relationships among H. formosa cytb haplotypes with nodal support using maximum-likelihood tree searches and bootstrap searches (500 pseudoreplicates) in GARLI 0.97 , partitioning the data by codon position, ((1 + 2), 3). We analyzed our multilocus dataset using similar GARLI runs partitioning the data into cytb codon-position subsets and an RPS7-gene subset, unlinking parameters across subsets. GARLI runs relied on DNA substitution models selected using the decision theory algorithm in DT-ModSel  (Supplement S1). We compared GARLI results to phylogenetic results from Bayesian coalescent-dating analyses (below) with similar site models. Because incomplete lineage sorting can obscure intraspecific phylogenetic relationships, we derived a 'species tree’ for our cytb haplotypes using the 'Minimize Deep Coalescences’ method . Although crude, this method increases probability of obtaining accurate population trees using even a single locus and yields insights into the degree of lineage sorting . We also inferred phylogenetic relationships among populations from the allozyme dataset. We imported unbiased Nei’s D estimates (above) into PAUP* 4.0b10  and used them to construct a midpoint-rooted neighbor-joining tree with total-distance branch lengths. We tested the genetic distinctiveness of the resulting clades in multivariate space using principal coordinates analysis (PCoA) conducted in GENALEX. Last, we inferred networks of cytb and RPS7 haplotype relationships with connections exceeding 95% parsimonious probability in TCS 1.21 .
We evaluated temporal diversification by estimating H. formosa-Poeciliidae sp. divergence, H. formosa population/clade divergence and other node ages using coalescent-dating analyses of our multilocus matrix in BEAST 1.74 . To ensure convergence, we ran three independent searches (MCMC = 108, sampled every 4000 generations; burn-in = 107) using identical priors and relaxed, uncorrelated lognormal molecular clocks. We linked tree models but partitioned the data into cytb codon-position subsets ((1 + 2), 3) and an RPS7-gene subset and unlinked clocks and parameters across subsets. Coalescent constant population size tree priors were used. We set uniform priors on rates drawn on by the lognormal relaxed clock that, for cytb, spanned protein-coding mtDNA gene substitution rates reported for teleost fishes ('fish rate’ = 0.017-0.14 × 10-8 subs/site/yr, per-lineage; refs. in [81, 82]), and for RPS7 spanned an arbitrary range of reasonable rates (1.0 × 10-10-1.0 × 10-8 subs/site/yr) consistent with nDNA rates for freshwater fishes (unpublished data). By including Poecilia (subgenus Limia) outgroups in our analysis, we were able to add a calibration point constraining the split between P. (L.) domicensis (from Cuba) and P. (L.) vittata (Hispaniola) to 17-14 Ma based on Poecilia (L.) phylogeny  and geological dates for the separation of Cuba and Hispaniola , following Doadrio et al. . We calibrated this node using a lognormal prior (mean in real space = 1, log standard deviation = 1.25, offset = 14). We used similar calibrations to model basal Pseudoxiphophorus divergence between 11-5 Ma, following Agorreta et al. , and to model the tree root age (Poeciliidae) to 39.9 Ma with an extended tail (log standard deviation = 2.0) based on the oldest fossils known for the family from the Paleocene-Eocene Maiz Gordo and Lumbrera formations, Argentina . We summarized posterior parameter distributions and calculated parameter-trace effective sample sizes (ESS) in TRACER (http://tree.bio.ed.ac.uk/software/tracer). We archived our sequence alignments and maximum-likelihood and Bayesian tree results in TreeBASE (Submission 14718; http://purl.org/phylo/treebase/phylows/study/TB2:S14718).
We used coalescent simulations of mtDNA in Mesquite 2.73  to statistically discriminate among phylogeographical scenarios representing the expansion-contraction model and two alternatives. This approach guards against the possibility that stochastic coalescences might make it more difficult to detect a signature of one or more models [43, 44]. Population tree models in the simulations (described in Results, Supplement S1) spanned a total time (tTotal) equal to tree depths of 3.741 × 106 generations, derived by converting the mean intraspecific time to the most recent common ancestor (tMRCA) estimate from BEAST from absolute time to coalescent time (=1.247 Ma/T). We scaled branching intervals so that they summed to tTotal. We estimated overall Ne for simulations (assumed equal to ancestral Ne) by summing female Ne (Nef) estimates from θ values, calculated for each SAMOVA-inferred population using the equation θ = 2Neμ (for mtDNA, μ = cytb rate estimated in BEAST). Thetas were derived from θW (above), and from Bayesian Θ estimates for SAMOVA groups, averaged from three replicate MIGRATE-N 3.1.3  runs each using 1 long chain (3 × 108) and adaptive heating (chains = 1, 1.5, 3, 104; Supplement S1). We also estimated effective numbers of female migrants per generation (Nefm) from mtDNA in MIGRATE-N. Populations (internal branches) were modeled using groupings of subpopulations (areas, SAMOVA groups) appropriate to each hypothesis, and we scaled population-lineage sizes (internal branch widths) to their respective proportions of overall Ne, so they summed to overall Ne at all time points . Population-lineage sizes were equally distributed among tip branches (subpopulations) evolving from them. We conducted hypotheses testing in Mesquite by simulating 1000 mtDNA gene genealogies within the expansion-contraction model under neutral coalescence. For these simulated gene trees, we calculated the 'number of deep coalescences’ (nDC), a measure of gene tree-population tree discordance arising from incomplete lineage sorting . We evaluated model fit by comparing nDC for the 'best’ cytb maximum-likelihood gene tree to that from the simulated gene trees (treated as rooted). We rejected the expansion-contraction model in favor of an alternative by a one-tailed significance test (α = 0.05 level) .
Ecological niche modeling
Highly predicted (>80%) present and LIG ranges overlapped in patches around Lake Pontchartrain, Florida’s Apalachicola Bay-Apalachee Bay and western peninsula (Hillsborough River), and the northwest Atlantic Coast range margin. Under more conservative suitability criteria (logistic thresholds), present-day and LIG-modeled ranges covered >90% of the modern range. There were virtually no false negatives, but we obtained minor false positives where H. formosa does not occur (e.g. above the Fall Line, in the Bahamas).
Statistical model-evaluation analyses and comparisons of variable contributions supported these conclusions. Based on AUC scores from 10 replicates, models performed well with high discriminatory ability (LGM minimum training presence logistic threshold = 0.079, AUC = 0.920 ± 0.016; LIG minimum training presence logistic threshold = 0.022, AUC = 0.915 ± 0.018). Average test model gain (>1.5) also indicated >4-fold better-than-random prediction. The principal environmental predictor of occurrence was precipitation during warmest quarter, followed by temperature seasonality (Additional file 1: Tables S3-S4). MESS diagrams showed near-coastal regions experienced more novel climates during the LGM than the LIG, relative to present day, causing model extrapolation (Figure 3A,B). Still, novel environments were not a major issue: predictor variables most responsible for environmental novelty (novel limiting environments), e.g. isothermality, contributed ≤2.1% to MAXENT model prediction.
Population structure and genetic diversity
We found substantial mean intraspecific mtDNA haplotype diversity (± standard deviation: Hd = 0.934 ± 0.006) but low mean nucleotide diversity (π = 0.0066 ± 0.009) and surprisingly shallow genetic divergence (mean cytb p-distance = 0.66%, range: 0-1.61%), all of which are consistent with a recent, rapid population expansion for the species. Using DnaSP, we calculated an empirical estimate of θW = 0.0087 for all H. formosa populations combined. SAMOVA-group estimates of population mutation rate parameter were moderate, with θW ranging 0.0031-0.0053 (Additional file 1: Table S5).
All three methods of spatial genetic analyses recovered three locations of genetic distinctions that defined four groups (K) within the species range (Figure 1). The groupings were strongly supported by bootstrap analysis (BARRIER BP = 55-99); significant spatially explicit among-group variance partitioning (SAMOVA P < 0.0001); and AMOVAs (Supplement S1). While SAMOVA and BARRIER analyses were not completely concordant, both algorithms recovered genetic barriers outside of the ecological niche modeling-inferred LGM refugium, and both inferred complex structuring between the Apalachicola delta and Aucilla River. These results are most consistent with leading-edge colonization or allele surfing (or other factors) during post-glacial range expansion. However, they do not firmly reject vicariance-northeast colonization: barriers were identified between WCP and ACP regions, but not FL and ACP; groups were divided north-to-south between the Suwannee and Hillsborough Rivers; and no ACP genetic barriers were inferred.
Most genetic barriers reflected breaks between drainages, and we also found significant among-drainage genetic partitioning in Florida (AMOVA: Among-group variation = 33.73%, ΦCT = 0.34, P < 0.0001). The most striking SAMOVA-inferred subdivision was isolation of this northern Florida area from all other populations (WCP, south Florida peninsula, ACP). However, within northern Florida, SAMOVA connected geographically disjunct lakes/ponds (collections 13-15, 17, and 20), Womack Creek (10), and Robinson Creek (30), and linked disjunct Newnan’s Lake (35) and Tate’s Hell sites (7-8). BARRIER recovered contiguous groups (Figure 1) and split the WCP from other areas (barrier ii); however, given disparities with our other DNA and allozyme results, and barrier ii’s placement in a sampling gap, we interpreted the WCP group as a manifestation of the limitations of our mtDNA sampling.
Mantel tests indicated that these conclusions were not unduly influenced by spatial autocorrelation in cytb variation , whether rangewide (r = -0.007, P = 0.492; regression R2 = 5 × 10-5, t = -0.11, P = 0.913) or only within Florida (r = -0.054, P = 0.391; regression R2 = 0.003, t = -0.74, P = 0.460). MtDNA isolation-by-distance was thus overall lacking, reflecting low diversity and alleles shared among regions/populations; for example, haplotype 22 was common to all three regions and most ACP populations shared haplotype 4 (Table 1 and S6). We noted a similar lack of nDNA isolation-by-distance, with RPS7 haplotype 2 present in each study region (Additional file 1: Table S7). These results agree with expansion-contraction predictions as well as the independently derived niche models.
Nonparametric tests of spatial-diversity predictions based on allozyme variation
Comparison (group 1-vs.-group 2)
Within-vs.-outside putative refuge
ACP-vs.-rest of range
WCP-vs.-rest of range
Within-vs.-outside putative refuge
WCP-vs. putative refuge
Neutrality and mismatch distribution tests for population expansions within Heterandria formosa
Expansion time (ka)
1.195 [0.000, 10.393]
24.5 [0.0, 213.4]
1.498 [0.000, 3.119]
30.8 [0.0, 64.0]
1.672 [0.000, 9.023]
34.3 [0.0, 185.3]
SAMOVA group 2
2.408 [1.256, 3.031]
49.4 [25.8, 62.2]
SAMOVA group 3
0.992 [0.492, 1.615]
20.4 [10.1, 33.2]
SAMOVA group 4
3.000 [0.410, 3.172]
61.6 [8.4, 65.1]
Phylogenetic relationships and coalescent-dating analyses
Several features of the sequence data supported predictions derived from the expansion-contraction model and were not consistent with the predictions of a vicariance-northeast colonization model. For one, no deep intraspecific east-west breaks were supported. This is not because there were no deep phylogenetic breaks detected; the predicted sister relationship between Poeciliidae sp. and H. formosa (maximum cytb p-distance = 9.01%) was strongly supported (Figure 6A, Additional file 1: Figure S3), indicating a deep split between Mexico and the Atlantic Plain. Our 'best’ cytb haplotype gene tree (ln L = -2000.8824) revealed shallow variation and a close Ochlockonee-Aucilla-Suwannee relationship, and supported differentiation (albeit incomplete) of SAMOVA groups 2 and 3. This topology differed greatly from the 'Minimize Deep Coalescences’ topology, indicating substantial incomplete lineage sorting (Figure 6B).
Allozyme phylogenetic analysis recovered four shallow, spatially overlapping clades that were also genetically distinct during PCoA multivariate projection, with axes 1 and 2 explaining 39.9% and 28.8% of the genetic variance, respectively (Figure 2). Consistent with expansion-contraction predictions, allozyme results corroborated SAMOVA-inferred population structure, with a north-south break between the Suwannee and Hillsborough Rivers, separating a group of northwestern plus south Florida populations from all others. Allozyme results also agreed with other phylogenetic results: the northwest Florida allozyme clade matched the strongly supported subclade 'a’ in Figure 6A. Contrasting our DNA sequence-based results, and consistent with vicariance-northeast colonization, neighbor-joining and PCoA recovered ACP populations as a genetically and geographically distinct clade/group. However, against vicariance-northeast colonization expectations, intraspecific clades/groups were surprisingly not split east-to-west near the Apalachicola River-region (Figure 2) or St. Mary’s River, as suggested by .
The model resulting from multilocus coalescent-dating in BEAST (Figure 6A; ln L ± standard error = -10,054.902 ± 0.118) yielded a geometric mean tMRCA of H. formosa samples dated to 1.247 Ma [0.620, 2.178] in the Early-Middle Pleistocene. The tMRCA estimate for subclade 'a’ was 0.377 Ma [0.153, 0.713]. Bayesian posterior distributions for these tMRCA estimates strongly rejected LGM divergence (P < 0.05; <5% of posterior samples estimated a tMRCA for this node younger than 22 ka). Thus, the deepest population structuring (interpopulation divergence) within H. formosa appears to have initiated before the LGM. We estimated divergence between Poeciliidae sp. and H. formosa at 8.478 Ma [3.939, 14.913] in the Late Miocene-Early Pliocene; confidence intervals for this split were much wider than the tMRCAs above but we still statistically rejected Pleistocene divergence (P < 0.05). Parameter-trace ESS scores >600 indicated thorough Markov chain mixing. Our BEAST runs estimated a geometric mean evolutionary cytb rate of 7.12 × 10-9 substitutions/site/yr, and an RPS7 rate of 1.01 × 10-9 substitutions/site/yr.
By summing regional group Nef estimates based on empirical mtDNA θW (Additional file 1: Table S1) and Θ (range: 0.00106-0.00360) estimates, we respectively inferred overall Nef ≈ 6.05 × 105 and 3.53 × 106. We also inferred migration rates among SAMOVA groups (mean Nefm range: 1.184-2.887); however, Bayesian posterior Nefm distributions peaked at the lower limit of resolution, with 95% confidence intervals including zero. Setting ancestral Ne equal to the Nef estimates above, we ran separate simulations testing three hypotheses modeled as hypothetical population trees (see Additional file 1: Supplement S1): (1) a null, 'fragmented ancestor’ model representing the expansion-contraction model, with a single size-constant ancestral population evolving through the Pleistocene before experiencing a 90% LGM reduction in Ne (range contraction), and subsequent colonization and diversification of local subpopulations (15-0 ka; summing to current Ne) while simultaneously expanding from a source population (south Florida refugium located by ecological niche modeling; Figure 3B); (2) a 'vicariance-northeast colonization’ model with a Early Pleistocene ancestor split basally east-west of Apalachicola River/Mobile Bay and ACP populations arising spontaneously post-LGM (15-0 ka) from a Florida refuge (colonized from a Saint Johns River source); and (3) a 'four-refugia’ vicariance model with SAMOVA population groups (Figure 1, results below) diverging from a Pleistocene ancestor. The number of deep coalescences for our 'best’ cytb gene tree fit into the population trees for each hypothesis was 9 for the null model, 14 for vicariance-northeast colonization, and 27 for four-refugia. Results of both simulation sets failed to reject the expansion-contraction hypothesis in favor of alternative models (P < 0.001), thus our gene tree is more consistent with a scenario of LGM population contraction and postglacial re-expansion than either of the alternative vicariance models.
Our results failed to reject most expansion-contraction model predictions (66.7% of tests; summarized in Additional file 1: Table S8). Indeed, several independent lines of evidence from geospatial and genetic data give concordant results consistent with this model. Here, we explore these results.
Ecological niche modeling suggests that H. formosa underwent range contraction to a single LGM refugium in the southwest Florida peninsula, an area that underwent less environmental-climatic change than the northern parts of the species range (Figure 3B). Paleoclimatic data, climate models and pollen records indicate southwest Florida maintained ~0-4°C lower temperatures, wetter tropical climates (+400-800 mm/yr annual precipitation) and mixed southern pine/hardwood forest to open woodland vegetation during the LGM [11, 27, 93], relative to today. By contrast, LGM conditions in northern Florida and the 'mainland’ Gulf-Atlantic plains were colder/drier with extreme ~8-10°C drops in winter temperatures, 400 mm/yr lower annual precipitation, and patchy northern forest and grasslands vegetation [27, 93, 94]. Combined with these data, our results strongly suggest H. formosa were extirpated from their northern range and survived the LGM in bioclimatically suitable peninsular areas, areas that, in fact, were the only parts of eastern North America that did not experience significantly colder winter and summer LGM temperatures .
While one paleobotanical synthesis suggested this putative refugial area was nothing more than arid sand-dune scrub at the time , other data  indicate that at this time south Florida was semiarid, containing suitable freshwater stream and wetland/marsh habitats where cold-intolerant freshwater organisms persisted. This deduction is consistent with pollen data showing aquatic macrophytes such as Brasenia, Isoëtes, Myriophyllum, Typha and Sparganium occurring within the Florida peninsula during the LGM .
Comparing the predicted LGM refuge of H. formosa with its much wider modern distribution (Figures 1 and 3C), we suggest that H. formosa rapidly recolonized most or all of its preceding range to the north (e.g. of the LIG, Figure 3A) during postglacial climatic 'amelioration.’ After 16-15 ka, eastern North American rivers underwent transitions from braided to meandering channel types, during which wetter-than-modern early Holocene conditions may have increased wetlands and river discharge . Conditions from this time onward may thus have favored expansion into parts of the species range, even though temperature and precipitation did not stabilize near their modern levels until somewhat later, ~9-6 ka . Such northward post-glacial expansions were unlikely to have been slowed by Younger Dryas cooling 12-11 ka (i.e. potentially negative effects on aquatic habitats), which did not affect areas below ~38°N latitude (, refs. therein). As a corollary, we reject vicariance-northeast colonization-predicted range dynamics and our results also seem to run counter to rangewide population stasis and connectivity expected under the southern crossroads hypothesis ; apparently, smaller areas of the North American Coastal Plain (south Florida) may have experienced prolonged ecological stability than in Mediterranean Basin coastal plains.
Ecological niche modeling assumes that relationships between present distributions, abiotic factors, and ecological interactions remain more or less at equilibrium, that species physiological limits have changed less than the magnitude of changes in abiotic variables, and that biotic features of a species’ niche affect historical distribution less than large-scale changes in abiotic conditions [18, 22, 23]. We cannot evaluate all of these assumptions, nor can we exclude possible bias due to coastline effects i.e. major influence of the Bermuda high-pressure zone on coastal precipitation, which was among the variables of highest (>30%) predictive importance (Additional file 1: Tables S3-S4). Nonetheless, our niche-based modeling results provide meaningful insights into H. formosa biogeography. Our paleoclimatic models yield highly accurate predictions, e.g. AUCs >0.90 and our results are consistent with existing ecophysiological data. For example, modern H. formosa populations reach critical thermal minima at ~10°C  and live close to this cold-season limit in northern parts of their range. Coinciding with this value, 88% of the H. formosa localities (Additional file 1: Figure S2) had mean temperature of coldest quarter values ≥10°C (±0.1°C), and this variable had sharply higher effects on MAXENT prediction above 10°C (JCB, unpublished data). Given cold-season temperatures across the species entire northern range declined to -8 to 4°C during the LGM [27, 93], it is reasonable to infer that glacial-stage extirpation occurred passively in these areas, with drought and colder winter temperatures reducing survival. In addition, the assumption of ecological niche conservatism on these spatial and temporal scales seems justified for this species. Differences in thermal tolerances among H. formosa populations from different parts of their range are minimal (~2 °C; ) compared to the temperature differences induced by climatic oscillations of the last glacial cycle above; thus differential physiological adaptation to temperature regimes likely has not had a confounding effect on our inferences. Our analyses of DNA sequence data and re-analyses of nuclear allozyme loci reinforce the paleoclimatic data to support expansion-contraction dynamics. In particular, the negative latitudinal cline suggests a glacial refugium in south Florida and northward recolonization [4, 9, 15]. Admixture in south Florida due to secondary contact of diverged lineages, i.e. following 'vicariance’ expected under vicariance-northeast colonization, or recent dispersal mediated by episodic connections (e.g. flooding, axial-valley connections) can be ruled out for two reasons. First, there is a distinct lack of phylogenetic divergence and, second, seven allozyme alleles in south Florida are private alleles not shared with populations in other areas. These data are also inconsistent with southern crossroads dynamics supported elsewhere, in which no clinal diversity patterns were witnessed . That our phylogeographical results also support expansion-contraction dynamics agreeing with the niche models also suggests that the assumption of ecological niche conservatism through time is essentially correct for H. formosa.
Of course, it is difficult to discount the possibility of 'microrefugia’  in small habitat pockets elsewhere. This is because our niche-based models predict responses to prevailing macroclimatic conditions, not microclimatic conditions below the data-layer grain (1 km2) or within local habitats. If any secondary refugium existed, the most likely candidate area would seem to have been the Tallahassee area, where we also found one private allozyme allele.
Several of our results agree with previous studies [42, 54]. Similar to Baer  and despite extreme distances of >2000 km between our collection localities and those in the allozyme dataset, we find H. formosa are characterized by shallow genetic divergences and a single population-lineage (e.g. mean mtDNA p-distance = 0.66%; Figures 6 and 7). Combined with high Hd but low π estimates, this suggests recent postglacial expansion with too limited time for recovery to allow accumulation of large sequence divergence , and/or high historical gene flow levels. Also congruent with , we infer no genetic barriers within the Atlantic Coastal Plain. This is consistent with Baer’s 'northeast colonization’ hypothesis (ACP populations colonized following an Early Pleistocene high-sea stand following 1.5 Ma), which he argued was also supported by lower allozyme diversity and a lack of ACP isolation-by-distance  (see below). Heterandria formosa also contains four significantly differentiated mtDNA population groups subdivided outside of the putative refugium, or between the refugium and recolonized areas (Figures 1, 2 and 7). The latter break, inferred between the Suwannee and Hillsborough Rivers (i.e. SAMOVA groups 1 vs. 2-4; Figures 1 and 7), corresponds roughly to a north-to-south break between Gulf-draining Waccasassa and Withlacoochee Rivers described by Baer . This break creates a clade of Louisiana plus south Florida samples in our phylogenetic and multivariate ordination analyses of the allozyme-frequency data, and Baer’s (Figure 2; though our results are based on different genetic distances). This Suwannee-Hillsborough break thus appears robust to different data types and methods and is shallow in both studies (e.g. ~0.6-0.8% mtDNA divergence) implying recently limited gene flow between these rivers. Perhaps more importantly, the Suwannee-Hillsborough break falls (with other genetic barriers) outside the predicted LGM refugium but is not correlated with deep phylogenetic divergence. This is consistent with the action of density-dependent processes, such as priority effects e.g. [2, 9], or allele surfing , operating during leading-edge recolonization. Given their patchy distribution [39, 42] and propensity for explosive population growth , it is reasonable to assume that H. formosa quickly fill new habitats once they arrive and that when local patches become extinct they are recolonized from nearby populations. Thus if northward postglacial expansion occurred through small numbers of colonists, alleles might have surfed along a relatively quickly expanding wave front. Our niche-based models suggest postglacial genetic divergence between these rivers may have partly been maintained by lack of predicted suitable habitat between these rivers (Figure 3).
By contrast, we find several points of nuclear-mitochondrial discordances, indicating differences between studies. For example, allozyme and mtDNA results yield mixed support for the expansion-contraction prediction that isolation-by-distance should occur within the refugium, but not outside of it, due to a longer demographic timeline and geographically restricted dispersal within the refuge. Allozymic variation exhibits nearly rangewide isolation-by-distance in our study and Baer’s . Baer  also inferred populations were at or near drift-migration equilibrium, using a 2-dimensional stepping stone model with cross-peninsula gene flow. Global Nem estimates were high (>4, range: 2.5-37.5), suggesting local H. formosa populations have not evolved independently but that gene flow has importantly shaped their evolution ; and in a second study, Nem within and across two peninsular Florida rivers ranged >3-400 . In Baer , ACP populations were exceptional, forming a relatively homogeneous monophyletic group with no isolation-by-distance consistent with recent colonization, consistent with our mtDNA and allozyme findings. However, contrasting these patterns, we find: a lack of mtDNA isolation-by-distance across the range; Nem suggesting zero on-going mtDNA gene flow (see Supplement S1); and ACP populations as indistinct from other regions (i.e. nonmonophyletic; Figures 2, 6 and 7). Incongruent with Baer’s  allozyme results, the historical signal of the mtDNA genome also does not indicate east-to-west genetic differentiation between WCP and ACP regions (Figure 1).
Such nuclear-mitochondrial discordances are traditionally explained by several factors. The first and most common explanation is balancing selection (e.g. during genetic hitchhiking) mainly on allozymes , some of which (particularly metabolic protein loci) are known to evolve non-neutrally. However, this is rejected outright by the data: Lewontin-Krakauer tests on the allozymes [42, 54] and multiple tests of mtDNA neutrality (e.g. Additional file 1: Table S5) reject the hypothesis of selection operating in either genome, as does the concordant lack of deep phylogenetic structure within either marker-class. A second explanation is that these discordances arise from the contrasting modes of transmission, evolution, and resolution of these marker-types; this could plausibly explain a lack of isolation-by-distance in the mtDNA, but isolation-by-distance in the allozymes. Compared with haploid, maternally transmitted mtDNA, allozymes are diploid and evolve roughly 10× slower with ~4-fold larger Ne and coalescence times (independent of mutation rate), and higher gene-flow rates [29, 101, 102]. As a consequence, mtDNA provide a higher-resolution view of population structure, e.g. at finer spatial scales (supported by our data), and more likely reflect homogenizing effects of gene flow hence lack of isolation-by-distance. By contrast, due to lower resolution, allozymes may display isolation-by-distance in northern areas even when expansion occurred from a single southern refugium, if postglacial expansion(s) did not involve leptokurtic (long-distance) dispersal e.g. reviewed in [2, 4, 9]. Regarding discordant migration patterns, allozyme Nem estimates and Ne estimates seem preferable in the case of H. formosa (mtDNA typically drastically underestimate Ne in high-gene flow species ). Still, their large coalescence times mean allozyme-inferred Nem may reflect averages across transmission routes over four times the history of mtDNA, rather than current processes. Thus, we interpret our results and Baer’s  as indicating that historical gene flow (high Nem) among populations nearing migration-drift equilibrium, and large Ne, has influenced H. formosa phylogeographical structuring. This could partly explain the low observed genetic divergences, but can we reconcile this scenario with the isolation-by-distance findings, taken at face value? Taking such a scenario as compatible with the expansion-contraction model assumes most population structure has arisen since the LGM and that sufficient time has passed for H. formosa allozyme loci to reach equilibrium. This is supported by the data if we assume Nm = 10 (similar to , or higher ) and that populations must have achieved migration-drift equilibrium since the LGM, or G = 57,000 generations (19 ka*T). Given our Nef (~6 × 105-3 × 106) and ~1.2 Ma intraspecific tMRCA estimates, the number of generations required to approach equilibrium is recent, ~60,829 generations, and since the LGM based on the equation G = 1/(2 m + 1/(2Ne)) . Thus, the above interpretation seems reasonable in light of population genetics theory; equilibrium could have arisen even faster in this species, as postglacial recolonization and population expansion can rapidly generate population structuring .
Another factor potentially explaining our results is sex-biased migration. Heterandria formosa exhibit strongly female-biased sex ratios due to higher male mortality rates, as demonstrated by predation trials  and otolith ring counts (JT, personal observation). In light of this, female dispersal/gene flow should outpace male dispersal even assuming equal dispersal propensities and rates between sexes. Thus, aside from differences in marker resolution, female-biased dispersal has likely also contributed importantly to nuclear-mitochondrial discordance above, specifically lack of mtDNA isolation-by-distance. This is because, under female-based dispersal, females carry both mtDNA and male/female nuclear genes out of populations while males are effectively philopatric (geographically restricted), thus we expect matrilineal alleles to more likely homogenize among localities while nuclear genomes remain site-dependent ; this would result in a relatively weaker matrilineal signal of isolation-by-distance. Unfortunately, we cannot reliably empirically evaluate effects of female-biased survival and dispersal on our other genetic results given comparable metrics influenced by geographical distance (e.g. linearized FST) are not estimable from the mtDNA and allozyme-frequency data, and low RPS7-gene sampling restricts possible mtDNA-nDNA comparisons. Last, consistent with the observation that deep coalescences appear to influence the H. formosa mtDNA/gene tree (e.g. Figure 6B), patterns of discordance between markers herein may to some extent reflect historical processes of incomplete lineage sorting in addition to the homogenizing effects of gene flow [significant when Nm > 4 (reviewed in ), as in H. formosa]. In future studies, sampling multiple, unlinked nuclear loci will be necessary to tease apart the relative contributions of these two opposing forces within H. formosa using coalescent-based models.
In addition to the geospatial and genetic analyses above, our phylogenetic and historical-demographic results also met predictions of the expansion-contraction model. While our inferences regarding ancestral derived/populations were limited by reduced south Florida sampling, phylogenetic and spatial patterns of ancestral/derived populations (alleles) between the phylogeny and network analyses congruently follow expectations derived from the model (e.g. Figure 4, Additional file 1: Table S8), and star-like patterns in the network suggest a history of population expansions. Our argument that H. formosa fits the expansion-contraction model is also reinforced by mismatch analyses and mtDNA genetic equilibrium tests (Fu’s FS, R2), which fail to reject demographic expansion models (e.g. for SAMOVA-inferred population groups; see text; Table 3, S5; Figure 5) and thus capture genetic signatures of predicted historical changes in Ne containing information about the record and timing of expansion [61, 70, 92]. The 95% confidence intervals on expansion parameter τ suggest the onset of expansions has overlapped the LGM in most population groups and regions analyzed, and lower 95% confidence intervals extending into the present may be indications that expansion is on-going or ceased only recently, e.g. in north Florida (Table 3; Figure 5). Because balancing selection seems not to have been at play (as per above), low-frequency peaks in our mismatch distributions likely reflect secondary Holocene expansions rather than selection (Figure 5).
These demographic results generally agree with genetic patterns in European taxa (despite possibly earlier onset of H. formosa expansions) that are inferred to have experienced postglacial population expansions e.g. [4, 17]; however, in our case, integrating niche-based modeling and genetic analyses provides an additional check that it is reasonable to conclude H. formosa expansion has been both spatial and demographic (Figure 3; Table 3). Interestingly, our mtDNA genetic equilibrium tests only point to recent population bottlenecking within WCP samples (Table 3) spanning extensive wetlands and swamps of the Mississippi River Delta, a major outlet of postglacial meltwater flow and alluvial deposition. This could reflect a small population surviving the LGM in one or more WCP microrefugia not identified in our paleoclimatic models; however, no evidence exists that a population-lineage experienced drift in isolation in the WCP since the LGM (e.g. Figure 6B). We thus interpret this bottleneck as a result of founder events during postglacial recolonization.
Although the distribution of H. formosa spans nearly the entire lowland Coastal Plain (Figure 1), its genetic patterns do not completely coincide with those described in several other species with overlapping present-day distributions. At least 20 species, including 50-67% of freshwater fishes and turtles examined to date, display pronounced east-west Apalachicola phylogeographical [29–32] and distributional  breaks. In four freshwater fishes—bowfin (Amia calva) and three sunfishes (Lepomis gulosus, Lepomis microlophus, and Lepomis punctatus)—this pattern has been linked to vicariant isolation in upper reaches of different drainage basins that were reduced in size by eustatic high-stands 50-80 m ASL of the Pliocene-Early Pleistocene interglacials [105, 106]. Thus, major phylogenetic divergence(s) should be expected in this region in as-yet unsampled freshwater taxa. Remarkably, however, neither patterns of multilocus phylogenetic structuring or mtDNA-population structuring, nor a coalescent simulation perspective on neutral mtDNA evolution, support such a pattern within H. formosa. Variation at protein loci in H. formosa also departs from this vicariance model, revealing shallowly diverged and geographically overlapping population groups/lineages without deep genetic breaks at predicted barriers. While the common pattern of Gulf drainages containing Atlantic-coast haplotypes (but not vice versa) is present in the allozymes (Figure 2), neither mtDNA nor allozyme data support clear east-west Apalachicola splits.
Our estimated tMRCA suggests intraspecific diversification of H. formosa has coincided with Pleistocene diversification within A. calva, which is codistributed with H. formosa in the Coastal Plain but ranges much more widely across eastern North America (including the Mississippi and St. Lawrence River basins). This indicates that recent isolation within constraints of Pleistocene-Holocene (rather than earlier) drainage geomorphology, e.g. drainage divides, has influenced genetic variation of both H. formosa (i.e. genetic distinctiveness of mtDNA at AMOVAs structured by Florida drainages) and A. calva. However, additional studies will be necessary to test whether these species diversification has actually been synchronous, and to evaluate finer-scale patterns relative to drainage basins. It would also be interesting to examine whether Coastal Plain freshwater fish species that are presently codistributed with H. formosa responded to Late Pleistocene disruptions in climate and sea-level with similar range-shifting responses, possibly involving refugia; however, the only species in which effects of Pleistocene expansion-contraction dynamics have been rigorously tested is H. formosa (this study). Here, the eastern/Atlantic intraspecific lineages (Florida peninsula-Atlantic seaboard) present in many freshwater taxa will likely present good candidates for testing for range expansion-contraction dynamics similar to H. formosa. For example, the eastern/Atlantic lineages recovered by Bermingham and Avise  exhibit <2-4% genetic divergences, indicating Pleistocene tMRCAs (assuming the standard 2%/Myr clock for vertebrate mtDNA ), and phylogroups or private alleles unique to the Florida peninsula. While ad hoc explanations that interspecific differences in the position of the Apalachicola break (particularly geographical distributions of the eastern lineages) resulted from unique histories of dispersal and range-shifting in response to Pleistocene climatic fluctuations have been advanced [29, 31, 32], the actual pattern of such dynamics remains an unanswered question. Additional studies similar to ours of previously studied fish taxa, e.g. in , as well as understudied aquatic plants, insects, freshwater fishes (e.g. Fundulus seminolis, Jordanella floridae, etc.), and crayfishes from the Florida peninsula will therefore permit testing the generality of the scenario of historical range-shift and population dynamics uncovered within H. formosa, as well as general predictions of the expansion contraction model , within the Gulf-Atlantic Coastal Plain freshwater biota.
Understanding historical responses of species to the climatic and sea-level fluctuations of the Late Pleistocene has emerged as a major goal of ecology, evolution and biogeography [2, 4, 15, 29]. We find the evolutionary history of Heterandria formosa indicates southward LGM range contraction and large-scale postglacial expansion, as predicted under the Pleistocene expansion-contraction model . Species-specific expansion-contraction dynamics may therefore have played a more general role in shaping the evolutionary history of Coastal Plain biota than previously thought. By adding to considerable variation surrounding common biogeographical themes recovered among even codistributed species in this region, our results bolster growing appreciation for the complexity of phylogeographical patterning in North America’s southern refugia. Our study demonstrates the benefits of taking an integrative approach to biogeographical hypotheses generation and testing (e.g. more realistic hypotheses testing through niche modeling and statistical phylogeography) and showcases how similar approaches can be used by biologists in the future to gain further insights into the evolutionary history of the Gulf-Atlantic Coastal Plain biota, as well as the important role of this region as a dynamic refuge during the last glacial cycle. Indeed, a better understanding of the tempo and mode of species responses to the Pleistocene glacial-interglacial cycles will likely aid predicting and managing ecosystem responses to present-day and future climate change .
Availability of supporting data
Sequences are deposited in GenBank (accession numbers: KF632895-KF633114, cytb, and KF633115-KF633133, RPS7). Sequence alignments and phylogenetic trees are available in the TreeBASE repository: http://purl.org/phylo/treebase/phylows/study/TB2:S14718.
We are grateful to M F. Breitman, E. Castro, and J. A. Wooten for helpful comments on this manuscript, and to A. L. Almendra for critical assistance with ecological niche modeling analyses. We thank P. J. Unmack for discussions and for graciously donating GIS data and samples from Hillsborough and St. Johns Rivers, Florida. Sampling was conducted under permits from Louisiana, Mississippi, Alabama, Georgia, Florida, and South Carolina issued to MS, JT and The University of Alabama Ichthyological Collection (UAIC). We thank A. Contreras-Balderas and G. Hubbard for assistance with field collections. Research was funded by Brigham Young University, including a Mentoring Environment Grant to JBJ that covered molecular sequencing costs; The University of Alabama; and The Florida State University.
- Gibbard PL, Head MJ, Walker MJC, and the Subcommission on Quaternary Stratigraphy: Formal ratification of the Quaternary System/Period and the Pleistocene Series/Epoch with a base at 2.58 Ma. J Quaternary Sci. 2009, 25: 96-102.Google Scholar
- Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405 (6789): 907-913. 10.1038/35016000.PubMedGoogle Scholar
- Hewitt GM: Ice ages: their impact on species distributions and evolution. Evolution on Planet Earth. Edited by: Rothschild LJ, Lister AM. 2003, London, UK: Academic Press, 339-361.Google Scholar
- Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philos T Roy Soc B. 2004, 359 (1442): 183-195. 10.1098/rstb.2003.1388.Google Scholar
- EPICA community members: Eight glacial cycles from an Antarctic ice core. Nature. 2004, 429: 623-628. 10.1038/nature02599.Google Scholar
- Gates DM: Climate Change and its Biological Consequences. 1993, Sunderland, MA: Sinauer AssociatesGoogle Scholar
- Lomolino M, Riddle B, Whittaker R, Brown J: Biogeography. 2010, Sunderland, MA: Sinauer Associates, 4Google Scholar
- Hof C, Levinsky I, Araújo MB, Rahbek C: Rethinking species’ ability to cope with rapid climate change. Glob Change Biol. 2011, 17: 2987-2990. 10.1111/j.1365-2486.2011.02418.x.Google Scholar
- Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996, 58 (3): 247-276.Google Scholar
- Bennett KD, Tzedakis PC, Willis KJ: Quaternary refugia of north European trees. J Biogeogr. 1991, 18 (1): 103-115. 10.2307/2845248.Google Scholar
- Jackson ST, Webb RS, Anderson KH, Overpeck JT, Webb T, Williams JW, Hansen BCS: Vegetation and environment in Eastern North America during the last glacial maximum. Quat Sci Rev. 2000, 19: 489-508. 10.1016/S0277-3791(99)00093-1.Google Scholar
- Davis MB: Quaternary history and the stability of forest communities. Forest Succession. Edited by: West DC, Shugart HH, Botkin DB. 1981, New York, NY: Springer-Verlag, 132-177.Google Scholar
- Delcourt PA, Delcourt HR: Vegetation maps for eastern North America: 40,000 yr. BP to the present. Geobotany II. Edited by: Romans RC. 1981, New York, NY: Plenum, 123-165.Google Scholar
- FAUNMAP Working Group: Spatial response of mammals to late quaternary environmental fluctuations. Science. 1996, 272 (5268): 1601-1606. 10.1126/science.272.5268.1601.Google Scholar
- Provan J, Bennett KD: Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008, 23 (10): 564-571. 10.1016/j.tree.2008.06.010.PubMedGoogle Scholar
- Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF: Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998, 7 (4): 453-464. 10.1046/j.1365-294x.1998.00289.x.PubMedGoogle Scholar
- Hewitt GM: Post-glacial recolonization of European biota. Biol J Linnean Soc. 1999, 68 (1–2): 87-112.Google Scholar
- Waltari E, Hijmans RJ, Peterson AT, Nyari AS, Perkins SL, Guralnick RP: Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions. PLoS One. 2007, 2 (7): e563-10.1371/journal.pone.0000563.PubMed CentralPubMedGoogle Scholar
- Gómez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography in Southern European Refugia: Evolutionary Perspectives on the Origins and Conservation of European Biodiversity. Edited by: Weiss S, Ferrand N. 2007, Dordrecht: Springer-Verlag, 155-188.Google Scholar
- Canestrelli D, Cimmaruta R, Nascetti G: Population genetic structure and diversity of the Apennine endemic stream frog, Rana italica—insights on the Pleistocene evolutionary history of the Italian peninsular biota. Mol Ecol. 2008, 17: 3856-3872. 10.1111/j.1365-294X.2008.03870.x.PubMedGoogle Scholar
- Porretta D, Canestrelli D, Urbanelli S, Bellini R, Schaffner F, Petric D, Nascetti G: Southern crossroads of the Western Palaearctic during the Late Pleistocene and their imprints on current patterns of genetic diversity: insights from the mosquito Aedes caspius. J Biogeogr. 2011, 38: 20-30. 10.1111/j.1365-2699.2010.02385.x.Google Scholar
- Svenning JC, Normand S, Kageyama M: Glacial refugia of temperate trees in Europe: insights from species distribution modelling. J Ecol. 2008, 96: 1117-1127. 10.1111/j.1365-2745.2008.01422.x.Google Scholar
- Fløjgaard C, Normand S, Skov F, Svenning JC: Ice age distributions of European small mammals: insights from species distribution modelling. J Biogeogr. 2009, 36: 1152-1163. 10.1111/j.1365-2699.2009.02089.x.Google Scholar
- Dépraz A, Cordellier M, Hausser J, Pfenninger M: Postglacial recolonization at a snail’s pace (Trochulus villosus): confronting competing refugia hypotheses using model selection. Mol Ecol. 2008, 17: 2449-2462. 10.1111/j.1365-294X.2008.03760.x.PubMedGoogle Scholar
- Otto-Bliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A, CAPE Last Interglacial Project members: Simulating arctic climate warmth and icefield retreat in the last interglaciation. Science. 2006, 311: 1751-1753. 10.1126/science.1120808.PubMedGoogle Scholar
- Chen JH, Curran HA, White B, Wasserburg GJ: Precise chronology of the last interglacial period: 234U-230Th data from fossil coral reefs in the Bahamas. Geol Soc Am Bull. 1991, 103: 82-97. 10.1130/0016-7606(1991)103<0082:PCOTLI>2.3.CO;2.Google Scholar
- Shin SI, Liu Z, Otto-Bliesner B, Brady EC, Kutzbach JE, Harrison SP: A simulation of the last glacial maximum climate using the NCAR-CCSM. Clim Dynam. 2003, 20 (2–3): 127-151.Google Scholar
- Lambeck K, Chappell J: Sea level change through the last glacial cycle. Science. 2001, 292: 679-686. 10.1126/science.1059549.PubMedGoogle Scholar
- Avise JC: Phylogeography: the History and Formation of Species. 2000, Cambridge, MA: Harvard University PressGoogle Scholar
- Soltis DE, Morris AB, McLachlan JS, Manos PS, Soltis PS: Comparative phylogeography of unglaciated eastern North America. Mol Ecol. 2006, 15 (14): 4261-4293. 10.1111/j.1365-294X.2006.03061.x.PubMedGoogle Scholar
- Avise JC: Molecular population structure and the biogeographic history of a regional fauna: a case history with lessons for conservation biology. Oikos. 1992, 63 (1): 62-76. 10.2307/3545516.Google Scholar
- Bermingham E, Avise JC: Molecular Zoogeography of Fresh-Water Fishes in the Southeastern United-States. Genetics. 1986, 113 (4): 939-965.PubMed CentralPubMedGoogle Scholar
- Swift CC, Gilbert CR, Bortone SA, Burgess GH, Yerger RW: Zoogeography of the fresh water fishes of the south-eastern United States: Savannah River to Lake Ponchartrain. Zoogeography of North American Freshwater Fishes. Edited by: Hocutt CH, Wiley EO. 1985, New York, NY: Wiley, 213-265.Google Scholar
- Ronquist F: Dispersal-vicariance analysis: a new approach to the quantification of historical biogeography. Syst Biol. 1997, 46 (1): 195-203. 10.1093/sysbio/46.1.195.Google Scholar
- Swenson NG, Howard DJ: Clustering of contact zones, hybrid zones, and phylogeographic breaks in North America. Am Nat. 2005, 166 (5): 581-591. 10.1086/491688.PubMedGoogle Scholar
- Edwards CE, Soltis DE, Soltis PS: Molecular phylogeny of Conradina and other scrub mints (Lamiaceae) from the southeastern USA: evidence for hybridization in Pleistocene refugia?. Syst Bot. 2006, 31: 191-205.Google Scholar
- Edwards CE, Soltis DE, Soltis PS: Using patterns of genetic structure based on microsatellite loci to test hypotheses of current hybridization, ancient hybridization and incomplete lineage sorting in Conradina (Lamiaceae). Mol Ecol. 2008, 17: 5157-5174. 10.1111/j.1365-294X.2008.03985.x.PubMedGoogle Scholar
- Parks CR, Wendel JF, Sewell MM, Qiu Y-L: The significance of allozyme variation and introgression in the Liriodendron tulipifera complex (Magnoliaceae). Am J Bot. 1994, 81: 878-889. 10.2307/2445769.Google Scholar
- Martin FD: Heterandria formosa Agassiz. Atlas of North American Freshwater Fishes. Edited by: Lee DS, Gilbert CR, Hocutt CH, Jenkins RE, McAllister DE, Stauffer JRJr. 1980, Raleigh, NC: North Carolina State Museum of Natural History, 547-Google Scholar
- Leips J, Travis J: The comparative expression of life-history traits and its relationship to the numerical dynamics of four populations of the least killifish. J Anim Ecol. 1999, 68 (3): 595-616. 10.1046/j.1365-2656.1999.00311.x.Google Scholar
- Schrader M, Travis J, Fuller RC: Do density-driven mating system differences explain reproductive incompatibilities between populations of a placental fish?. Mol Ecol. 2011, 20: 4140-4151. 10.1111/j.1365-294X.2011.05264.x.PubMedGoogle Scholar
- Baer CE: Species-wide population structure in a southeastern U.S. freshwater fish, Heterandria formosa: gene flow and biogeography. Evolution. 1998, 52 (1): 183-193. 10.2307/2410933.Google Scholar
- Knowles LL, Maddison WP: Statistical phylogeography. Mol Ecol. 2002, 11 (12): 2623-2635. 10.1046/j.1365-294X.2002.01637.x.PubMedGoogle Scholar
- Nielsen R, Beaumont MA: Statistical inferences in phylogeography. Mol Ecol. 2009, 18 (6): 1034-1047. 10.1111/j.1365-294X.2008.04059.x.PubMedGoogle Scholar
- Knowles LL: Statistical phylogeography. Annu Rev Ecol Evol Syst. 2009, 40: 593-612. 10.1146/annurev.ecolsys.38.091206.095702.Google Scholar
- Carstens BC, Richards CL: Integrating coalescent and ecological niche modeling in comparative phylogeography. Evolution. 2007, 61 (6): 1439-1454. 10.1111/j.1558-5646.2007.00117.x.PubMedGoogle Scholar
- Richards CL, Carstens BC, Knowles LL: Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. J Biogeogr. 2007, 34 (11): 1833-1845. 10.1111/j.1365-2699.2007.01814.x.Google Scholar
- Pluzhnikov A, Donnelly P: Optimal sequencing strategies for surveying molecular genetic diversity. Genetics. 1996, 144: 1247-1262.PubMed CentralPubMedGoogle Scholar
- Felsenstein J: Accuracy of coalescent likelihood estimates: do we need more sites, more sequences, or more loci?. Mol Biol Evol. 2006, 23: 691-700.PubMedGoogle Scholar
- Hrbek T, Seckinger J, Meyer A: A phylogenetic and biogeographic perspective on the evolution of poeciliid fishes. Mol Phylogenet Evol. 2007, 43 (3): 986-998. 10.1016/j.ympev.2006.06.009.PubMedGoogle Scholar
- Unmack PJ, Bagley JC, Adams M, Hammer MP, Johnson JB: Molecular phylogeny and phylogeography of the Australian freshwater fish genus Galaxiella, with an emphasis on dwarf galaxias (G. pusilla). PLoS One. 2012, 7 (6): e38433-10.1371/journal.pone.0038433.PubMed CentralPubMedGoogle Scholar
- Chow S, Takeyama H: Intron length variation observed in the creatine kinase and ribosomal protein genes of the swordfish Xiphias gladius. Fish Sci. 1998, 64: 397-402.Google Scholar
- Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25 (11): 1451-1452. 10.1093/bioinformatics/btp187.PubMedGoogle Scholar
- Baer CF: Population structure in a south-eastern US freshwater fish, Heterandria formosa. II. Gene flow and biogeography within the St. Johns River drainage. Heredity. 1998, 81: 404-411. 10.1046/j.1365-2540.1998.00410.x.Google Scholar
- Phillips SJ, Anderson RP, Schapire RE: Maximum entropy modeling of species geographic distributions. Ecol Model. 2006, 190 (3–4): 231-259.Google Scholar
- Elith J, Graham CH, Anderson RP, Dudik M, Ferrier S, Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, Lehmann A, Li J, Lohmann LG, Loiselle BA, Manion G, Moritz C, Nakamura M, Nakazawa Y, Overton JMC, Peterson AT, Phillips SJ, Richardson KS, Scachetti-Pereira R, Schapire RE, Soberón J, Williams S, Wisz MS, Zimmerman NE: Novel methods improve prediction of species' distributions from occurrence data. Ecography. 2006, 29 (2): 129-151. 10.1111/j.2006.0906-7590.04596.x.Google Scholar
- Stockwell DRB, Peterson AT: Effects of sample size on accuracy of species distribution models. Ecol Model. 2002, 148 (1): 1-13. 10.1016/S0304-3800(01)00388-X.Google Scholar
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25 (15): 1965-1978. 10.1002/joc.1276.Google Scholar
- Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA, Chang P, Doney SC, Hack JJ, Henderson TB, Kiehl JT, Large WG, McKenna DS, Santer BD, Smith RD: The Community Climate System Model version 3 (CCSM3). J Clim. 2006, 19: 2122-2143. 10.1175/JCLI3761.1.Google Scholar
- Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 28 (10): 2731-2739. 10.1093/molbev/msr121.PubMed CentralPubMedGoogle Scholar
- Excoffier L, Foll M, Petit RJ: Genetic consequences of range expansions. Annu Rev Ecol Evol Syst. 2009, 40: 481-501. 10.1146/annurev.ecolsys.39.110707.173414.Google Scholar
- Monmonier MS: Maximum-difference barriers: an alternative numerical regionalization method. Geogr Anal. 1973, 3: 245-261.Google Scholar
- Manni FE, Guerard E, Heyer E: Geographical patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by "Monmonier's algorithm". Hum Biol. 2004, 76: 173-190. 10.1353/hub.2004.0034.PubMedGoogle Scholar
- Excoffier L, Lischer HEL: Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010, 10 (3): 564-567. 10.1111/j.1755-0998.2010.02847.x.PubMedGoogle Scholar
- Hood GM: PopTools, Version 3.0.3. 2012, Available at http://www.cse.csiro.au/poptoolsGoogle Scholar
- Dupanloup I, Schneider S, Excoffier L: A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002, 11 (12): 2571-2581. 10.1046/j.1365-294X.2002.01650.x.PubMedGoogle Scholar
- Peakall R, Smouse PE: GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006, 6: 288-295. 10.1111/j.1471-8286.2005.01155.x.Google Scholar
- Mantel N: The detection of disease clustering and a generalized regression approach. Cancer Res. 1967, 27: 209-220.PubMedGoogle Scholar
- Wakeley J: The effects of subdivision on the genetic divergence of populations and species. Evolution. 2000, 54 (4): 1092-1101.PubMedGoogle Scholar
- Harpending HC: Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994, 66 (4): 591-600.PubMedGoogle Scholar
- Hammer Ø, Harper DAT, Ryan PD: PAST: paleontological statistics software package for education and data analysis. Palaeontol Electron. 2001, 4 (1): 1-9.Google Scholar
- Mcdonald JH, Kreitman M: Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991, 351 (6328): 652-654. 10.1038/351652a0.PubMedGoogle Scholar
- Fay JC, Wu CI: Hitchhiking under positive Darwinian selection. Genetics. 2000, 155 (3): 1405-1413.PubMed CentralPubMedGoogle Scholar
- Zwickl DJ: Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. PhD thesis. 2006, Austin: The University of Texas at Austin,Google Scholar
- Minin V, Abdo Z, Joyce P, Sullivan J: Performance-based selection of likelihood models for phylogeny estimation. Syst Biol. 2003, 52 (5): 674-683. 10.1080/10635150390235494.PubMedGoogle Scholar
- Maddison WP, Knowles LL: Inferring phylogeny despite incomplete lineage sorting. Syst Biol. 2006, 55 (1): 21-30. 10.1080/10635150500354928.PubMedGoogle Scholar
- Knowles LL, Carstens BC: Estimating a geographically explicit model of population divergence. Evolution. 2007, 61 (3): 477-493. 10.1111/j.1558-5646.2007.00043.x.PubMedGoogle Scholar
- Swofford D: PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods), version 4.0. 2001, Sunderland, MA: Sinauer AssociatesGoogle Scholar
- Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9 (10): 1657-1659. 10.1046/j.1365-294x.2000.01020.x.PubMedGoogle Scholar
- Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.PubMed CentralPubMedGoogle Scholar
- Waters JM, Burridge CP: Extreme intraspecific mitochondrial DNA sequence divergence in Galaxias maculatus (Osteichthys: Galaxiidae), one of the world's most widespread freshwater fish. Mol Phylogenet Evol. 1999, 11 (1): 1-12. 10.1006/mpev.1998.0554.PubMedGoogle Scholar
- Burridge CP, Craw D, Fletcher D, Waters JM: Geological dates and molecular rates: fish DNA sheds light on time dependency. Mol Biol Evol. 2008, 25 (4): 624-633. 10.1093/molbev/msm271.PubMedGoogle Scholar
- Hamilton A: Phylogeny of Limia (Teleostei: Poeciliidae) based on NADH dehydrogenase subunit 2 sequences. Mol Phylogenet Evol. 2001, 19 (2): 277-289. 10.1006/mpev.2000.0919.PubMedGoogle Scholar
- Iturralde-Vinent MA: Meso-Cenozoic Caribbean paleogeography: implications for the historical biogeography of the region. Int Geol Rev. 2006, 48: 791-827. 10.2747/0020-6818.104.22.1681.Google Scholar
- Doadrio I, Perea S, Alcaraz L, Hernandez N: Molecular phylogeny and biogeography of the Cuban genus Girardinus Poey, 1854 and relationships within the tribe Girardinini (Actinopterygii, Poeciliidae). Mol Phylogenet Evol. 2009, 50: 16-30. 10.1016/j.ympev.2008.09.014.PubMedGoogle Scholar
- Agorreta A, Dominguez-Dominguez O, Reina RG, Miranda R, Bermingham E, Doadrio I: Phylogenetic relationships and biogeography of Pseudoxiphophorus (Teleostei: Poeciliidae) based on mitochondrial and nuclear genes. Mol Phylogenet Evol. 2013, 66: 80-90. 10.1016/j.ympev.2012.09.010.PubMedGoogle Scholar
- Pascual R, Bond M, Vucetich MG: El Subgrupo Santa Bárbara (Grupo Salta) y sus vertebra- dos. Cronología, paleoambientes y paleobiogeografía. 1981, San Luis, Argentina: VIII Congreso Geologico Argentino, Actas III, 746-758.Google Scholar
- Maddison WP, Maddison DR: Mesquite: a modular system for evolutionary analysis. Version 2.73. 2010, Available from http://mesquiteproject.orgGoogle Scholar
- Beerli P, Felsenstein J: Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. P Natl Acad Sci USA. 2001, 98 (8): 4563-4568. 10.1073/pnas.081068098.Google Scholar
- Shepard DB, Burbrink FT: Phylogeographic and demographic effects of Pleistocene climatic fluctuations in a montane salamander. Plethodon fourchensis. Mol Ecol. 2009, 18 (10): 2243-2262. 10.1111/j.1365-294X.2009.04164.x.Google Scholar
- Maddison WP: Gene trees in species trees. Syst Biol. 1997, 46 (3): 523-536. 10.1093/sysbio/46.3.523.Google Scholar
- Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9 (3): 552-569.PubMedGoogle Scholar
- Prentice LC, Bartlein PJ, Webb T: Vegetation and climate change in eastern North America since the last glacial maximum. Ecology. 1991, 72: 2038-2056. 10.2307/1941558.Google Scholar
- Watts WA: Late-Quaternary vegetation history at White Pond on the Inner Coastal-Plain of South-Carolina. Quat Res. 1980, 13: 187-199. 10.1016/0033-5894(80)90028-9.Google Scholar
- Sawada M, Viau AE, Gajewski K: The biogeography of aquatic macrophytes in North America since the last glacial maximum. J Biogeogr. 2003, 30: 999-1017. 10.1046/j.1365-2699.2003.00866.x.Google Scholar
- Leigh DS: Terminal Pleistocene braided to meandering transition in rivers of the Southeastern USA. Catena. 2006, 66: 155-160. 10.1016/j.catena.2005.11.008.Google Scholar
- Baer CF, Travis J: Direct and correlated responses to artificial selection on acute thermal stress tolerance in a livebearing fish. Evolution. 2000, 54 (1): 238-244.PubMedGoogle Scholar
- Rull V: Microrefugia. J Biogeogr. 2009, 36: 481-484. 10.1111/j.1365-2699.2008.02023.x.Google Scholar
- Leips J, Travis J, Rodd FH: Genetic influences on experimental population dynamics of the least killifish. Ecol Monogr. 2000, 70: 289-309. 10.1890/0012-9615(2000)070[0289:GIOEPD]2.0.CO;2.Google Scholar
- Karl SA, Avise JC: Balancing selection at allozyme loci in oysters: implications for nuclear RFLPs. Science. 1992, 256: 100-102. 10.1126/science.1348870.PubMedGoogle Scholar
- Brown WM: Evolution of animal mitochondrial DNA. Evolution of Genes and Proteins. Edited by: Nei M, Koehn RK. 1983, Sunderland, MA: Sinauer Associates, 62-88.Google Scholar
- Moore WS: Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear-gene trees. Evolution. 1995, 51: 627-629.Google Scholar
- Crow JF, Aoki K: Group selelction for a polygenic behavioral trait: estimating the degree of population subdivision. Proc Natl Acad Sci USA. 1984, 81: 6073-6077. 10.1073/pnas.81.19.6073.PubMed CentralPubMedGoogle Scholar
- Richardson JML, Gunzburger MS, Travis J: Variation in predation pressure as a mechanism underlying differences in numerical abundance between populations of the poeciliid fish Heterandria formosa. Oecologia. 2006, 147 (4): 596-605. 10.1007/s00442-005-0306-y.PubMedGoogle Scholar
- Cronin TM, Szabo BJ, Ager TA, Hazel JE, Owens JP: Quaternary climates and sea levels of the U.S. Atlantic Coastal Plain. Science. 1981, 211 (4479): 233-240. 10.1126/science.211.4479.233.PubMedGoogle Scholar
- Riggs SR: Paleoceanographic model of Neogene phosphorite deposition, US Atlantic continental margin. Science. 1983, 223: 123-131.Google Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.