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 N
m 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, N
m 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; N
m 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 N
e 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 N
m estimates and N
e estimates seem preferable in the case of H. formosa (mtDNA typically drastically underestimate N
e in high-gene flow species ). Still, their large coalescence times mean allozyme-inferred N
m 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 N
m) among populations nearing migration-drift equilibrium, and large N
e, 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 N
ef (~6 × 105-3 × 106) and ~1.2 Ma intraspecific t
MRCA 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/(2N
e)) . 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 F
ST) 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 F
2), 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 N
e 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 t
MRCA 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 t
MRCAs (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.