Skip to main content

Extremophile Poeciliidae: multivariate insights into the complexity of speciation along replicated ecological gradients



Replicate population pairs that diverge in response to similar selective regimes allow for an investigation of (a) whether phenotypic traits diverge in a similar and predictable fashion, (b) whether there is gradual variation in phenotypic divergence reflecting variation in the strength of natural selection among populations, (c) whether the extent of this divergence is correlated between multiple character suites (i.e., concerted evolution), and (d) whether gradual variation in phenotypic divergence predicts the degree of reproductive isolation, pointing towards a role for adaptation as a driver of (ecological) speciation. Here, we use poeciliid fishes of the genera Gambusia and Poecilia that have repeatedly evolved extremophile lineages able to tolerate high and sustained levels of toxic hydrogen sulfide (H2S) to answer these questions.


We investigated evolutionary divergence in response to H2S in Gambusia spp. (and to a lesser extent Poecilia spp.) using a multivariate approach considering the interplay of life history, body shape, and population genetics (nuclear miscrosatellites to infer population genetic differentiation as a proxy for reproductive isolation). We uncovered both shared and unique patterns of evolution: most extremophile Gambusia predictably evolved larger heads and offspring size, matching a priori predictions for adaptation to sulfidic waters, while variation in adult life histories was idiosyncratic. When investigating patterns for both genera (Gambusia and Poecilia), we found that divergence in offspring-related life histories and body shape were positively correlated across populations, but evidence for individual-level associations between the two character suites was limited, suggesting that genetic linkage, developmental interdependencies, or pleiotropic effects do not explain patterns of concerted evolution. We further found that phenotypic divergence was positively correlated with both environmental H2S-concentration and neutral genetic differentiation (a proxy for gene flow).


Our results suggest that higher toxicity exerts stronger selection, and that divergent selection appears to constrain gene flow, supporting a scenario of ecological speciation. Nonetheless, progress toward ecological speciation was variable, partially reflecting variation in the strength of divergent selection, highlighting the complexity of selective regimes even in natural systems that are seemingly governed by a single, strong selective agent.


Darwin’s [1] view that speciation can be thought of as a continuum has gained increasing traction in recent years (e.g., [210]). Within this conceptual framework, pairs of populations move bidirectionally along a continuum between panmixis and complete reproductive isolation: populations can either progress towards speciation or collapse back into panmixis even after substantial population divergence (e.g., [1115]). Several factors shape the speciation continuum and determine where a given population pair falls along it. All else equal (i.e., assuming that the same speciation mechanisms act in replicate population pairs), progress towards complete speciation should increase as a function of time, with populations undergoing divergent selection or genetic drift for longer time periods moving farther along the continuum toward total reproductive isolation than populations that have diverged more recently. However, ‘all else’ is usually not equal, and so differences in progress towards complete speciation can also be explained by additional or stronger selective agents at play in only some replicate populations, mutation-order effects, differences in the genetic architecture of the founder populations, and/or stochastic effects (e.g., [7, 8, 16]).

The comparative approach, using multiple populations across replicated environmental gradients, not only provides a tool for uncovering mechanisms underlying speciation in natural systems—i.e., how reproductive isolation arises—but also patterns of convergent and non-convergent trait evolution (e.g., [7, 8, 1719]). Nonetheless, even comparative studies usually do not trace progression toward speciation in a single lineage through time but rather use comparisons among closely related taxa that have reached different stages of divergence, the underlying assumption being that a single lineage would follow a similar course through time as speciation unfolds (see also discussion in [20]). Studies employing such a comparative approach have usually focused on single traits or character suites, exemplified by several studies on variation in only life history or morphology (e.g., [2124]). This neglects the possibility that multiple character suites may evolve in concert (e.g., [2527]). In fact, only by simultaneously investigating multiple types of traits in conjunction with pertinent environmental features can we more fully understand the complexity of adaptive diversification and how speciation unfolds in natural systems [7, 8, 25, 2729]. In our present study, we use livebearing fishes (family Poeciliidae) of the genera Gambusia (mosquitofishes) and Poecilia (mollies) to assess divergence in multiple phenotypic character suites (i.e., life history and morphology) across a geographically replicated ecological gradient created by toxic hydrogen sulfide (H2S) (Fig. 1 and Table 1). We further measure the degree of neutral genetic differentiation between geographically adjacent population pairs as a proxy for gene flow—and thus, as an estimate of where population pairs reside along the speciation continuum.

Fig. 1
figure 1

a Overview of the general study area. (B-E) Detailed views of collection sites for (b) Gambusia affinis in Oklahoma, (c) Gambusia holbrooki in Florida, (d) Gambusia hubbsi from northern Andros Island, The Bahamas, and (e) Gambusia eurystoma and Gambusia sexradiata in Tabasco, southern Mexico. Numbers correspond to sampling sites as described in Table 1. Yellow: sulfidic habitats; blue: non-sulfidic sites. Panel (a) was created using GOOGLE EARTH VS. (©2015 Google Inc., Mountain View, CA, USA), while panels (be) were created using GOOGLE MAPS (Map data ©2015 Google, INEGI)

Table 1 Overview of sample sites and sample sizes used for analyses

To make specific predictions about the evolution of life histories and morphology, we need a clear understanding of the various sources of selection in sulfidic habitats. H2S toxicity results predominantly from its interference with mitochondrial bioenergetics, which inhibits oxidative phosphorylation [30, 31]. The suppression of aerobic respiration is aggravated by the reactivity of H2S leading to extreme hypoxia in aquatic environments; nevertheless, a number of livebearing fishes (Poeciliidae) have successfully colonized naturally sulfidic habitats [31]. Adaptation to H2S has been studied especially within the genus Poecilia from southern Mexico, where evolutionarily independent lineages have colonized sulfidic springs in four river drainages (e.g., [9, 3240]). Living in H2S-rich environments should be accompanied by life-history trait modifications, particularly through selection on offspring size [41]. First, as the ratio of sulfide influx to sulfide oxidation is crucial for efficient detoxification, larger neonates are expected to exhibit higher tolerance to elevated environmental H2S concentrations due to a lower body surface-to-volume ratio [42, 43]. Second, larger neonates are expected to have a competitive advantage in sulfidic habitats, which are thought to be resource-limited [31, 44, 45]. Congruent with these considerations, extremophile poeciliids of nine different species are characterized by larger offspring size at birth, and as a result of the classic life-history trade-off between offspring size and fecundity [46], are further characterized by a reduced fecundity compared to their closest relatives from non-sulfidic waters [41]. The question remains, however, whether larger offspring size in extremophile poeciliids is indeed a ubiquitous response to exposure to even low concentrations of environmental H2S, or whether extremophile populations could be found that lack this adaptation.

With respect to morphology, divergence in southern Mexico is mainly the result of extremophile Poecilia having larger heads coupled with an increased gill surface area to enhance respiratory capacities—and thus survival—in hypoxic H2S-springs [32, 34, 47]. Besides its direct impact on survival, relative head size may also provide a cue mediating mate choice, as female Poecilia in southern Mexican streams prefer males from their own population over immigrant (sulfide-adapted) males [9, 33].

Here, we start by focusing on all mosquitofish populations known to inhabit sulfidic waters. Specifically, we sampled four different clades of mosquitofish (G. affinis, G. holbrooki, G. hubbsi, and the two sister species G. sexradiata and G. eurystoma) to address the question of how predictable and convergent phenotypic divergence is in response to H2S among Gambusia species (question 1). Based on physiological studies, life-history theory, and previous work on extremophile Poecilia spp. (see above), we predicted extremophile mosquitofish to have larger but fewer offspring, and to develop larger heads compared to their closest relatives from non-sulfidic waters. We also examined other traits (e.g., adult fat content, reproductive allocation, overall body shape) that might diverge as a result of resource limitation (e.g., less body fat in sulfidic habitats) and higher population densities (e.g., more streamlined body shapes in sulfidic habitats), but for which previous studies on extremophile poeciliids from southern Mexico did not always find convergent and predictable patterns [24].

We then integrated the mosquitofish data with data collected on extremophile Poecilia spp. [9, 24, 40, 47, 48] to address two additional questions: First, is the strength of any observed shift in one of the character suites (life history or morphology) correlated with divergence in the other (question 2)? If developmental interdependencies, genetic linkage, or pleiotropy play important roles for phenotypic diversification in extremophile Gambusia, we would expect correlations between these character suites both within (i.e., on the inter-individual level) and between populations. Alternatively, the presence of H2S could create multifarious selective regimes that act on suites of traits not necessarily developmentally or genetically integrated. In that case we would predict trait correlations between body shape variation and important life-history traits (e.g., offspring size and fecundity) between, but not within, populations.

Moreover, we asked whether divergent selection between toxic and benign environments generally results in reduced gene flow (estimated through a population genetic approach using nuclear microsatellites), and if stronger phenotypic divergence is typically correlated with greater reductions in gene flow (question 3)? For the first time, we thus approach ecological speciation in H2S-rich environments from a multivariate perspective that considers the interplay of life history, morphology, and population genetics. We focused on both shared and unique patterns of population divergence (or the lack thereof) and provide a unifying framework within which variable progress along the speciation continuum can be interpreted.


Life-history divergence

Adult life histories

In the ANOVA on SL using the subset of data from Gambusia spp., the factors sex, clade, and site(clade × H2S), as well as the interactions of sex × clade, sex × H2S, and sex × clade × H2S had a significant influence, while the presence of H2S itself was not significant (Table 2a). Based on η p 2, the most important predictors for size differences were sex, clade and site(clade × H2S) (Table 2a). Males were generally smaller than females (estimated marginal means, EMM ± s.e., males: 20.85 ± 0.17 mm; females: 27.49 ± 0.15 mm; sex-effect in Table 2a); however, clades also differed in size (clade-effect) and within each clade, populations from the same habitat type also differed strongly from each other (effect of site(clade × H2S)).

Table 2 Results of nested univariate analysis of variance (ANOVA) and multivariate analyses of covariance (MANCOVA) examining life-history and body shape variation of Gambusia spp. from four different clades, comparing sulfidic springs versus non-sulfidic habitats. The site(clade × H2S)-term was designated a random effect in all models. F-ratios in the MANCOVAs were approximated using Wilks’ values, partial variance was estimated using Wilks’ partial η 2, and terms with a partial variance at least half as strong as the strongest effect are given in bold for each model. Asterisks indicate terms that were tested using the mixed-model approach (see text)

In the mixed-model nested MANCOVA on adult life histories, the covariate SL, the factors sex, H2S, clade, and site(clade × H2S), as well as the interactions of SL × clade, sex × clade, sex × H2S, clade × H2S, and sex × clade × H2S had a significant influence on the examined traits (Table 2b). Based on our measure of effect size (η p 2), SL and sex had by far the strongest influence on adult life histories while the factor H2S was of minor importance, and no consistent pattern of divergence due to the presence of H2S was uncovered (Table 2b; Fig. 2a). All three life-history traits generally increased with increasing SL (Pearson correlation, lean weight: r p = 0.94, P < 0.0001; fat content: r p = 0.18, P < 0.0001; GSI: r p = 0.65, P < 0.0001). Compared to similar-sized females, males weighed more (EMM ± s.e. at SL = 25.19 mm, lean weight, males: 0.083 ± 0.001 g, females: 0.067 ± 0.001 g), but had less body fat (males: 5.23 ± 0.44 %, females: 6.45 ± 0.32 %), and had a smaller relative investment into reproduction (GSI, males: 1.68 ± 0.37 %, RA, females: 17.70 ± 0.27 %; Table 2b and Additional file 1: Table S3). Lean weight, fat content, and GSI responded inconsistently to the presence of H2S across clades and sexes (Additional file 1: Table S3).

Fig. 2
figure 2

a Adult life-history, (b) offspring-related life-history, and (c) body shape divergence between Gambusia spp. populations from sulfidic habitats (yellow) and non-sulfidic sites (blue). Depicted are mean divergence vector scores (± s.e.; derived from the H2S term in the MANCOVAs) for each site; these describe the linear combination of dependent variables exhibiting the greatest difference between sulfidic and non-sulfidic habitats in Euclidean space (please see Additional file 1: OSM 3 for details). Numbers correspond to sites as described in Table 1

Descriptive statistics (site-specific means) for all N = 730 individuals analyzed for life history traits are summarized in Additional file 1: Tables S1 and S2.

Offspring-related life histories

In the mixed-model nested MANCOVA on offspring-related life histories, the covariates SL and stage of development, the factors H2S, clade, and site(clade × H2S), as well as the interactions clade × H2S, SL × clade, SL × H2S, stage of development × clade, and stage of development × H2S had a significant influence on the evaluated traits (Table 2c); however, based on η p 2, the presence of H2S and SL had by far the strongest influence (Table 2c). As predicted, the H2S-effect was mainly due to differences in embryo lean weight and fecundity (Fig. 2b), with females in sulfidic habitats producing significantly larger but fewer offspring (EMM ± s.e. for SL = 28.07 mm and stage of development = 27.87, embryo lean weight: non-sulfidic = 1.52 ± 0.04 mg, sulfidic = 3.04 ± 0.03 mg; fecundity: non-sulfidic = 14.98 ± 0.75 mg, sulfidic = 9.09 ± 0.66 mg; Additional file 1: Table S3). Fecundity increased with increasing female SL, but embryo fat content and embryo lean weight did not vary with female size (fecundity: r p = +0.67, P < 0.0001; embryo fat content: r p = -0.06, P = 0.24; embryo lean weight: r p = -0.08, P = 0.10).

Morphological divergence

In total, we examined N = 914 individuals of Gambusia spp. from 17 sites (Table 1). Centroid size, sex, clade, H2S, site(clade × H2S), and the interaction of clade × H2S had a significant influence on body shape (Table 2d). According to η p 2, the strongest differences were uncovered between clades and sexes, but body shape also differed between fish from sulfidic and non-sulfidic habitats (Table 2d). Thin-plate spline transformation grids visualizing the shape differences along the divergence vector derived from the H2S-term indicate that extremophile mosquitofish tended to have larger heads, more slender bodies, elongated caudal peduncles, and more posteriorly positioned dorsal fins (Fig. 2c). The degree of population differentiation varied among clades, with the strongest divergence found in G. affinis and the G. eurystoma/sexradiata complex (clade-effect), while G. hubbsi exhibited a trend toward slightly smaller heads and deeper bodies in the sulfidic habitats (clade × H2S-effect; Fig. 2c). The primary difference between sexes was in the position of the anal fin (sex-effect), which is modified into a copulatory organ (gonopodium) in males and positioned more anteriorly than the female anal fin [49].

Population genetic differentiation

Descriptive statistics for site-specific means of genetic variability of the N = 383 genotyped individuals are provided in Additional file 1: Table S5. We found variable degrees of genetic differentiation between fish from sulfidic and non-sulfidic habitats in G. holbrooki, ranging from virtual panmixis (F ST = 0.006) to moderate genetic structure (F ST = 0.249; Fig. 3a-c). In the Green Springs system, we detected K = 2 as the uppermost hierarchical level of population structure, with individuals from population 5 (Green Springs) being assigned to one genetic cluster and individuals from both non-sulfidic populations 6 and 7 (Lake Monroe and Flood Drain) being assigned to another cluster (Fig. 3a). However, the situation in the Florida Panhandle is considerably different. Using all sampling sites in one analysis we detected K = 2 clusters with animals originating from population 9 (Ditch Highway 98) being assigned to one and all other individuals (regardless of the presence of H2S) being assigned to the other cluster (Additional file 1: Figure S2). To infer the effect of H2S on gene flow between adjacent habitats, we decided to split our Florida Panhandle data and conducted separate, independent analyses. For the westernmost cluster of sites, populations 8 and 9 (Panacea Mineral Springs and Ditch Highway), we detected K = 3 to be the most likely number of clusters, involving multiple clusters within one site (Fig. 3b). A test for migration or hybridization within the last two generations was conducted using the USEPOPINFO model implemented in STRUCTURE, but could only detect minor probability for recent migration events between both populations (Q < 0.09). For the cluster of populations 10 to 12 (Newport Springs, St. Marks River, and Ditch St. Marks), the most likely level of population structure was obtained for K = 1, thus pointing towards unimpeded gene flow between sulfidic and non-sulfidic sites (Fig. 3c).

Fig. 3
figure 3

Population assignment using STRUCTURE version 2.3.2. Sulfidic and surrounding non-sulfidic habitats of (ac) G. holbrooki, (d) G. eurystoma/sexradiata, and (e) G. hubbsi. K = 1 was recovered as the most likely number of genetic clusters in c, K = 2 in a and e, and K = 3 in b and d; given are individual relative assignment values in yellow (sulfidic habitats; left side) as well as blue and olive (non-sulfidic sites; right side), divided by a black line and sorted by the relative assignment score for each population separately. Numbers correspond to sites as described in Table 1

K = 3 was uncovered as the uppermost hierarchical level of population structure in Mexico. The sulfidic site inhabited by G. eurystoma (population 13) formed a genetic cluster distinct from all G. sexradiata populations (populations 17–21; Fig. 3d). Pairwise F ST-values supported the results obtained from STRUCTURE, being highest in all pairwise comparisons of G. eurystoma and G. sexradiata (F ST ≥ 0.227).

The separation between sulfidic and non-sulfidic sites was also clear for the Bahamas sites, where K = 2 was uncovered as the uppermost hierarchical level of population structure: G. hubbsi from sulfidic Archie’s Blue Hole (population 23) formed one genetic cluster that was distinct from the nearby non-sulfidic habitats (populations 24 and 26; Fig. 3e).

The partial Mantel test on the global dataset uncovered a trend (albeit non-significant) for stronger differentiation between divergent habitats than between similar habitats (r P = 0.11, one-tailed P = 0.067), while controlling for differences among clades.

Joint evolution of life histories, morphologies, and reproductive isolation?

Life-history divergence in offspring-related traits was correlated with morphological divergence across all extremophile poeciliid complexes (N = 10, r S = 0.58, one-tailed P = 0.041; Fig. 4a). We further found an individual-level correlation between offspring-related life histories and body shape for G. holbrooki females from one of Florida’s three sulfide springs (population 10, Newport Springs: N = 44, r P = 0.45, one-tailed P = 0.001; Fig. 4b), but a similar correlation was lacking from Florida’s other two sulfide-spring populations (population 5, Green Springs: N = 33, r P = -0.10, one-tailed P = 0.29; population 8, Panacea Mineral Springs: N = 26, r P = -0.11, one-tailed P = 0.30, Fig. 4c-d; all three populations combined: N = 103, r P = -0.003, one-tailed P = 0.49). Congruently, a post-hoc ANCOVA on life-history divergence vector scores that included body shape divergence vector scores as a covariate revealed that the slopes for life-history divergence vs. morphological divergence differed significantly between the three populations as indicated by the significant interaction effect (site: F 2,97 = 120.15, P < 0.0001; body shape: F 1,97 = 0.31, P = 0.58; site × body shape: F 2,97 = 3.64, P = 0.03).

Fig. 4
figure 4

Scatterplots of body shape vs. life-history divergence vector scores (derived from the H2S-term in the MANCOVAs) for (a) each sulfidic complex, and (bd) pregnant female G. holbrooki from the three sulfide springs (sites 5, 8 and 10) in Florida. (e) Relationship between pairwise F ST (as a proxy for gene-flow, see Additional file 1: OSM 6) and phenotypic divergence (data for G. affinis not available) and (f) relationship between phenotypic divergence and H2S-concentrations (data for G. hubbsi not available). Black circles represent sulfidic systems inhabited by Gambusia, grey circles represent sulfidic systems inhabited by Poecilia; numbers correspond to sites as described in Table 1. Regression lines and 95 % CIs in grey are shown for significant correlations (see main text)

We also uncovered a significant correlation between phenotypic divergence and our proxy for gene flow (population genetic differentiation, F ST: N = 9, r S = -0.63, one-tailed P = 0.035; Fig. 4e), and the degree of phenotypic divergence for each sulfidic site increased with increasing H2S-concentration (N = 9, r S = 0.85, one-tailed P = 0.002; Fig. 4f).


Our present study leveraged replicated ecological gradients inhabited by poeciliid fishes of the genera Gambusia and Poecilia, and examined the speciation continuum using a multivariate approach by simultaneously considering life-history, morphological, and genetic differentiation, and the degree of reproductive isolation in different populations that have colonized H2S-rich waters. We uncovered strong signals of both shared and unique aspects of divergence in response to toxic H2S across lineages and populations. Specifically, we found (1) some broad patterns of predictable and convergent evolution of offspring-related life histories and adult body shape, (2) potential concerted evolution of different character suites, and (3) that reproductive isolation might emerge as a consequence of local adaptation to sulfidic environments. In other words, this study found evidence consistent with adaptation constraining gene flow, as trait divergence was negatively correlated with gene flow and positively correlated with environmental H2S concentrations, and genetic divergence increased between adjacent sulfidic and non-sulfidic populations in 4 of the 5 Gambusia systems examined here. To the extent that increased genetic divergence represents reduced gene flow (see Additional file 1: OSM 6 and Figure S6)—a scenario supported by the geographical proximity of the populations compared here—these results suggest that stronger divergent selection drives reduced gene flow in sulfidic systems inhabited by poeciliid fishes (see also [9]). However, we cannot fully exclude the alternative, but not mutually exclusive, direction of causation, whereby gene flow is partially responsible for constraining phenotypic divergence (e.g., [50, 51]).

We did not find any evidence for convergence across clades and toxic environments in adult life histories. Unlike offspring life histories and body shape, adult life-history traits do not seem to experience consistent directional selection from H2S, or at least do not exhibit consistent responses to such selection [24, 38]. Instead, the significant effects of sex, clade and site(clade × H2S) suggest that variation in adult life histories is likely shaped by other, often population-specific, factors, such as resource availability, sexual conflict and sexual selection or demography, or it reflects genetic drift [52].

Shared and unique aspects of phenotypic divergence and neutral genetic differentiation in extremophile Gambusia

We predicted extremophile mosquitofish to have larger but fewer offspring, and to develop larger heads compared to their closest relatives from non-sulfidic habitats. We uncovered strong overall convergence of morphologies and offspring-related life histories. Across four clades and six sulfur systems, we generally found phenotypic divergence consistent with our a priori predictions of local adaptation to H2S-rich environments (question 1).

There were, however, notable exceptions (see Additional file 1: OSM 6 for discussion of some minor patterns of unique divergence). Specifically, extremophile G. holbrooki from the Florida Panhandle (populations 8–12) displayed predictable, albeit weak, morphological divergence, but little to no life-history and neutral genetic divergence [i.e., population means for offspring size at birth and fecundity did not differ between populations from sulfidic and non-sulfidic waters in this region (Additional file 1: Table S1)]. The significant effect of H2S on offspring-related life histories in our multivariate models, however, indicates that at earlier developmental stages offspring size does indeed differ, but different degrees of embryonic weight loss result in more or less equally sized neonates. This is in contrast to a recent study demonstrating convergence of these traits in sulfidic waters in different lineages of poeciliid fishes [41]. Several non-mutually exclusive mechanisms could account for this observation: First, gene flow in this system might constrain local adaptation, as the spring complex at Panacea Mineral Springs (population 8) is only about 4 × 4 m wide before it drains into a large river where H2S immediately gets diluted; this habitat might be too small to allow for a self-sustaining, locally-adapted population (support for this comes from the fact that in the summer of 2015 we only caught P. latipinna here but no G. holbrooki). Second, H2S concentrations at Panacea Mineral Springs were relatively low (see Additional file 1: OSM 1) and might not be sufficient to drive local adaptation and the evolution of reproductive isolation. This hypothesis would also be congruent with the observed positive correlation between phenotypic divergence and H2S-concentration, as well as the negative correlation between phenotypic divergence and rates of gene flow. Nonetheless, this hypothesis alone does not explain why a similar pattern was observed at Newport Springs (population 10), which has H2S-concentrations apparently high enough to facilitate strong population divergence in G. affinis from Oklahoma (populations 1–4) as well as Poecilia mexicana from southern Mexico (Additional file 1: OSM 1; Fig. 4f) [9, 24, 34]. Third, while knowledge of the initial stages of ecological speciation in sulfidic habitats remains scarce, one possible scenario outlines a stepwise process: Initial colonizers might move back and forth between sulfidic and non-sulfidic waters, and gradually, subsequent generations of colonizers spend more and more time in sulfidic waters until populations become permanent residents. If this scenario were true, then it is likely that morphological changes evolve first because they are known to have a direct impact on survival of the colonizing individual, while offspring could still be birthed outside of the toxic waters during the early stages of colonization. Given the shallow genetic structure in these particular sulfur systems (see below), one could indeed argue that populations in the Florida Panhandle (populations 8–12) represent the very early stages of diversification in response to H2S. Fourth, speciation reversal (i.e., movement towards panmixis) can be initiated when environmental gradients break down or are overridden by another, stronger selective force (e.g., [1113, 15]). The Florida Panhandle is routinely impacted by tropical storms/hurricanes [53], and accompanying large-scale flooding may lead to temporary population admixture that resets population divergence at regular intervals. However, we consider this explanation to be least likely because the East Coast of Florida and the Bahamas are equally affected by these storms [54], yet population divergence is much stronger. Additionally, we found that patterns of genetic differentiation in sulfidic systems in southern Mexico were not affected by a catastrophic flooding event [36].

Concerted evolution of life histories and morphologies in Gambusia and Poecilia

We found that the magnitude of divergence in life history and morphology co-varied with one another across the ten sulfidic complexes (i.e., between-population correlation), and on an individual level both trait suites also co-varied in one out of three G. holbrooki sulfide spring populations (i.e., within-population correlation). Furthermore, overall phenotypic divergence also increased with increasing neutral genetic differentiation (a proxy for gene flow) as well as with increasing environmental H2S-concentrations (questions 2 and 3).

Positive within-population correlations between life-history and morphological divergence should be indicative of genetic linkage, pleiotropic effects, or developmental interdependencies as driving forces behind the apparent concerted evolution of both character suites; however, we only uncovered this correlation in one out of three populations from sulfidic waters for which data for both character suites were available for the same individuals. While this suggests that genetic linkage, pleiotropy or developmental interdependencies do not play a pivotal role in driving the observed patterns in body shape and offspring-related life histories, future studies using more population replicates will have to investigate this further.

Given the between-population correlation between morphological and life-history divergence, we argue that the presence of H2S creates—and is correlated with—complex environmental gradients. These exert multifarious selection that acts strongly and predictably on both genetically- and developmentally-independent character suites (see also [55]). We found additional support for this interpretation, because mean phenotypic divergence across study sites increased as a function of average H2S-concentrations.

Why exactly does the strength of phenotypic divergence correlate so strongly with H2S-concentration? First, our data indicate that population divergence in sulfidic habitats does not involve major threshold effects (with the possible exception of the comparatively low H2S-concentrations present at Panacea Mineral Springs, population 8), but rather mirror a gradual increase in divergent selection with increasing H2S-concentration. Our results leave open the possibility that a linear increase in phenotypic divergence occurs with increasing H2S concentration until a sufficiently high concentration is reached, after which phenotypes (and presumably selection) change little (i.e., possible asymptotic relationship in Fig. 4e). However, this is difficult to assess in nature, as sites with intermediate H2S concentrations are not currently known. Second, the observed phenotypic changes in toxic habitats likely come with associated costs. For example, females experience selection to maximize offspring number per clutch [56], which, all else being equal (i.e., barring other morphological changes that increase body cavity space), is constrained by selection to increase offspring size in sulfidic sites at least to the minimum size needed for efficient H2S-detoxification at the given H2S-concentration [4143]. On the morphological side, larger heads are known to reduce fast-start locomotor performance in poeciliid fishes [57, 58], which reduces escape ability and increases vulnerability to fast-striking predators like birds and fishes, and bird predation is known to be strong in many sulfidic habitats [59]. Thus, the morphological changes known to improve respiration in toxic habitats simultaneously render these sulfide-adapted fishes more vulnerable to predators.

Concerted evolution of phenotypic differentiation and reproductive isolation?

We found that population pairs with greater phenotypic differences exhibited reduced evidence of gene flow (estimated via genetic differentiation), and that phenotypic differences increased with increasing environmental H2S-concentrations. Therefore, when focusing on the interplay between phenotypic trait divergence, gene flow, and environmental H2S-concentrations, our data suggests that progress along the speciation continuum in extremophile poeciliids is generally promoted by increasing strength of divergent selection associated with higher H2S concentrations. In other words, higher H2S concentrations lead to stronger and more predictable phenotypic divergence, which in turn results in increased population genetic differentiation (and decreased gene flow), and thus shifts populations further along the speciation continuum towards complete reproductive isolation. This is in line with previous findings in extremophile Poecilia suggesting that local adaptations in life-history traits and morphology contribute to reproductive isolation, because they ought to result in a higher relative fitness of locally adapted resident fish compared to immigrants of the alternative ecotype (after [18]; e.g., [9, 33, 60]). This would also support our hypothesis (see above) that selection at Panacea Mineral Springs might simply be too weak to elicit pronounced population divergence. Nonetheless, future experimental data on pre- and post-mating isolation would likely help strengthen these results and yield new insights into the types of reproductive isolation promoting genetic divergence in each system.

For species and populations inhabiting environments characterized by H2S-concentrations ranging between ca. 30 and 40 μM (i.e., G. affinis from Vendome Well and G. holbrooki from Green Springs and Newport Springs in Florida; populations 1, 5, and 10, respectively), we uncovered substantial variation in the degree of phenotypic divergence and neutral genetic differentiation, indicating population-specific variation in the progress towards complete reproductive isolation. Several non-mutually exclusive scenarios are plausible that might help explain these differences. For example, if all of these populations do in fact experience the same speciation mechanisms, then populations undergoing divergent selection or genetic drift for longer should have moved further along the speciation continuum than populations that diverged more recently [7, 8, 16]. Furthermore, while it is tempting to assume that divergence in H2S-rich waters takes place in parapatry (based on the present situation found in sampled habitats), we cannot exclude the possibility that most systems that exhibit advanced stages of phenotypic divergence and neutral genetic differentiation have undergone at least temporary phases of allopatry during which populations from surrounding non-sulfidic waters temporarily collapsed. Evidence for this comes from southern Mexican sulphur mollies (Poecilia sulphuraria), which are genetically more closely related to northern Mexican P. mexicana limantouri than the southern Mexican P. mexicana mexicana they currently share the same continuous habitat with [34, 40, 61]. Moreover, H2S-discharge might be more temporally variable in some populations than others, with advanced divergence only being present in habitats that experience more constantly elevated H2S-concentrations. Finally, additional selective agents may be at play in some populations that either counteract or enhance H2S-induced selection, or this could result from mutation-order speciation, differences in the genetic architecture of the founder populations, and/or stochastic effects [7, 8, 16].


Our study highlights the diversity and complexity of organismal responses to immediate and evolutionary exposures to stressors, despite the fact that some physicochemical stressors (like H2S) appear to have clear-cut and predictable physical, chemical, or physiological effects on biological and non-biological processes at all levels of organization. This suggests an environmental complexity associated with the presence of physicochemical stressors, which can cause multifarious selective regimes through direct and indirect effects on abiotic and biotic components of an ecosystem. Such environmental complexity interacting with phenotypic integration and variation in genomic architecture can result in a pattern of shared and unique organismal responses similar to the one revealed here for extremophile poeciliids (with emphasis on Gambusia spp.) inhabiting sulfidic waters. Hence, integrative analyses of environmental and phenotypic complexity promise to gain a new level of sophistication in our understanding of life in extreme environments.


Study sites and sample collections

We collected mosquitofishes and mollies in sulfidic and adjacent non-sulfidic habitats in Oklahoma (G. affinis), Florida (G. holbrooki), Mexico (G. sexradiata, G. eurystoma, P. mexicana, and P. sulphuraria), and The Bahamas (G. hubbsi; Fig. 1; Table 1; see Additional file 1: Online Supplementary Material (OSM) 1 and 2 for details). All specimens designated for life-history and morphological analyses were preserved in 10 % formaldehyde solution, while specimens/fin clips collected for molecular genetic analyses were preserved in 70–95 % ethanol. For four species, we collected populations from both sulfidic and non-sulfidic habitats; however, this was not possible for G. eurystoma and P. sulphuraria, which are endemic to sulfide spring complexes in Tabasco, southern Mexico [62]. Both were therefore compared to their closely related sister species from non-sulfidic waters surrounding the habitats of the sulfide-spring endemics: G. eurystoma was compared to G. sexradiata [41, 63], and P. sulphuraria was compared to P. mexicana [40, 41]. Within each system, the degree of physical (geographical) separation between divergent habitats (i.e., sulfidic and non-sulfidic) is similar to, or frequently less than, that between similar (i.e., non-sulfidic) habitats (Fig. 4). Therefore, geographic isolation is not confounded with habitat type. Morphological data for G. affinis, P. mexicana, and P. sulphuraria were extracted from photographs used for previously published studies [40, 47]; however, we changed the original landmark configuration to match the one applied to all Gambusia spp. in the present study (Additional file 1: Figure S1). Offspring size at birth and fecundity data for most Gambusia and Poecilia populations were re-analyzed from [24, 41] (see Table 1 for details). All other data were collected specifically for this study.

Life-history and morphometric analyses

Dissections to collect male, female, and offspring-related life-history traits followed well-established protocols (e.g., [37, 38, 64]. We collected the following male and female life-history traits: standard length (SL [mm]), dry weight (g), lean weight (g), fat content (%), and reproductive investment [%; for males: testis dry weight divided by the sum of reproductive tissue dry weight and somatic dry weight (gonadosomatic index, GSI); for females: offspring dry weight divided by the sum of offspring dry weight plus somatic dry weight (reproductive allocation, RA)]. We note that GSI for males and RA for females do not capture total investment into reproduction, as it obviously ignores other costs of reproduction, such as energetic costs related to searching for mates, sneak-mating, and intrasexual competition for mates [65]. This is particularly true for males, where the relative behavioral costs of reproduction may exceed those captured by GSI [66]. Thus, patterns in total reproductive costs may not necessarily mirror those estimated via GSI and RA.

For females we also collected data on fecundity (number of developing offspring), offspring dry weight (mg), offspring lean weight (mg), and offspring fat content (%). Prior to statistical analyses we log10-transformed (male and female SL, male and female lean weight, and embryo dry and lean weight), square root-transformed (fecundity), or arcsine(square root)-transformed (male and female fat content, male and female GSI, embryo fat content) all life-history variables, and conducted subsequent z-transformation to meet assumptions of statistical analyses (i.e., these transformations facilitated normality of model residuals).

To quantify morphological differentiation, lateral photographs were taken using Canon DSLR cameras with a macro lens, and we then digitized 13 landmarks on each image (Additional file 1: Figure S1) using the software program tpsDig2 [67]. We performed geometric morphometric analysis as described in [32, 34] (see Additional file 1: OSM 3 for details), resulting in seven principal component axes (= relative warps) explaining 91.5 % of the morphological variance as variables for the statistical analyses.

Phenotypic differentiation based on life histories and body shape

We first tested for differences in adult body size (SL) between populations by conducting mixed-model analysis of variance (ANOVA) that included the following independent variables: clade (four levels: G. affinis, G. holbrooki, G. eurystoma/G. sexradiata, and G. hubbsi), sex, H2S (present or absent), and “site nested within clade-by-H2S” [random effect, hereafter: site(clade × H2S)]. In this and all subsequent statistical models we initially included all potential two-way and three-way interactions, but removed terms from the final model in a step-wise process if P > 0.1, with the exception that a term with P > 0.1 would be retained if a higher-order interaction term involving that term had a P < 0.1. We then conducted three separate mixed-model multivariate analyses of covariance (MANCOVA) examining variation in adult life history, offspring life history, and body morphology (see below for model structure). Assumptions of multivariate normal error distribution and homogeneity of variances and covariances were met for all analyses performed. Statistical significance was determined using an F-approximation from Wilks’s lambda for all model terms with the exception that we used an F-test using restricted maximum likelihood and the Kenward-Roger degrees of freedom adjustment [68] for clade, H2S, and clade × H2S to appropriately test these fixed effects while treating site as a random term (i.e., effectively treating site as the unit of replication; see also [69, 70]). The latter significance test was conducted using the MIXED procedure in SAS (SAS Institute, Cary, NC; sample code in the appendix of [69]), while all other tests were conducted in JMP (SAS). To quantify the relative importance of model terms, we calculated effect size using Wilks’s partial eta squared (η p 2) and calculated the relative variance as the partial variance for a given term divided by the maximum partial variance value in that model.

The first model tested for differentiation based on adult life histories and included lean weight, fat content, and GSI as dependent variables. To control for multivariate allometry, standard length (SL) was added as a covariate, and we included clade, sex, and H2S as fixed factors, and site(clade × H2S) as a random effect. The second model tested for differentiation based on offspring-related life histories and included fecundity, embryo lean weight, and embryo fat content as dependent variables. The fixed factors for this model were clade and H2S, while site(clade × H2S) was included as a random effect, and the covariates were SL and ‘embryonic stage of development’ (see [39] for details on embryo stages). The third model tested for phenotypic differentiation among sites based on body shape variation, using the seven retained relative warps as dependent variables. We tested for effects of centroid size to control for multivariate allometry and included clade, sex, and H2S as fixed factors, and site(clade × H2S) as a random effect. Running the analyses of adult life histories and body shape variation for both sexes separately confirmed our results and interpretations presented here (results not shown).

To assess the nature and strength of convergent life-history and morphological divergence in response to H2S among clades and sites, we performed a canonical analysis of the H2S-term of each MANCOVA to derive divergence vectors following [71]. Each divergence vector describes the linear combination of dependent variables that exhibits the greatest difference between habitat types in Euclidean space, while controlling for other terms in the model (see Additional file 1: OSM 4 for details). We visualized shape variation described by the H2S term included in the MANCOVA with thin-plate spline transformation grids using tpsRegr [72].

Population genetic analyses

We used 11 nuclear microsatellite loci to genotype N = 382 fish from 17 sites (Table 1; Fig. 4) in all species except G. affinis (for which no alcohol-preserved material was available, see Additional file 1: OSM 1) using previously established protocols [73, 74] (see Additional file 1: OSM 5 for details). We used ARLEQUIN v 3.5 [75] to calculate pairwise F ST-values between populations in each drainage and to calculate standard indicators of genetic variability (see Additional file 1: OSM 5). STRUCTURE v 2.3.4 [76] was employed to identify the number of genetically distinct clusters (K) in each drainage with the method outlined by [77] using the web-based tool STRUCTURE HARVESTER v 0.6.8 [78]. Ten iterations per K were run using the admixture model with a burn-in period of 106 generations, followed by 106 iterations for K = 1 up to twice the number of sampling sites included per area. Each simulation was performed using an ancestry model incorporating admixture, a model of correlated allele frequencies, and no prior information on locations.

To test whether divergent selection between sulfidic and non-sulfidic sites was associated with population genetic structure (i.e., higher F ST) while controlling for phylogenetic differences between clades, we employed a partial Mantel test with pairwise F ST-values obtained from FSTAT v 2.9.3 [79] as dependent variable, habitat difference as independent variable (0 = same habitat type, 1 = different habitat type), and clade as the covariate (0 = same clade, 1 = different clade). For this analysis, we compiled a global dataset including all Gambusia species for which microsatellite data were available (i.e., without G. affinis). To standardize comparisons across clades, we only included up to a maximum of two non-sulfidic sites for each sulfidic system in the first model, so that the final dataset consisted of the populations 5–12 for G. holbrooki, 13 and 16–17 for G. eurystoma/G. sexradiata, and 23–24 and 26 for G. hubbsi (see Table 1 for more details). Due to the lack of neutral genetic differentiation in the populations from the Florida Panhandle (see results), we conducted a second partial Mantel test that excluded those populations.

Joint evolution of life histories, morphologies, and reproductive isolation?

To increase sample size and statistical power, most of the following analyses were conducted on fishes from both genera (Gambusia and Poecilia; Fig. 4), thus adding offspring-related life-history data (N = 215), body shape data (N = 539), and population genetic data (N = 174) for Poecilia spp. from four additional sulfide spring systems (Table 1; Additional file 1: OSM 2 and Table S6 for details).

To investigate patterns of correlated phenotypic divergence across sulfide spring systems, we calculated average offspring-related life-history and morphological divergence scores for each sulfide spring complex (i.e., average divergence vector score for the non-sulfidic habitats – divergence vector score for the sulfide spring) and compared them with a Spearman rank correlation. We only included these two character suites as our results indicated that they exhibited the strongest shared responses to H2S (see results section). For Poecilia spp., we calculated scores along the divergence vectors derived for Gambusia spp. (described above) by projecting those individuals onto the relevant multivariate axes using the eigenvector coefficients.

We further investigated if the correlation between morphological and offspring-related life-history divergence can also be found on the individual level. For this we used a subset of the data, for which we had collected both life-history and morphological data from the same individuals [i.e., this included only pregnant G. holbrooki females from the three sulfide springs in Florida (sites 5, 8 and 10)]. The predicted positive correlation between variables (here and in subsequent correlations) was tested using one-tailed tests; non-parametric tests we used where appropriate to account for non-normally distributed data.

To create a single variable of phenotypic divergence, we subjected average life-history and morphological divergence scores for each sulfide system to a PCA based on a correlation matrix, and one principal component with an eigenvalue of 1.88 was retained, accounting for 94.0 % of the total variance. To investigate the relationship between overall phenotypic divergence and reproductive isolation, we calculated a Spearman rank correlation between these PC1 scores and pairwise F ST-values between a given sulfide spring and the geographically closest non-sulfidic site (F ST-values represented a good proxy for gene-flow in our system; see Additional file 1: OSM 6 for details). Due to a lack of population genetic data, G. affinis was excluded from this analysis.

Finally, we asked if the degree of ecologically-based divergent selection (approximated by average, system-specific H2S-concentrations) predicts the degree of phenotypic divergence. Since no data on H2S-concentrations for the G. hubbsi complex from the Bahamas were available, that clade was excluded from the analysis. We then tested for a correlation between PC scores for total phenotypic divergence (see above) and average H2S-concentrations of that particular sulfidic site using another Spearman rank correlation.


  1. Darwin C. On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life. London: John Murray; 1859.

    Book  Google Scholar 

  2. Mallet J, Beltrán M, Neukirchen W, Linares M. Natural hybridization in heliconiine butterflies: the species boundary as a continuum. BMC Evol Biol. 2007;7:28.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Hendry AP, Bolnick DI, Berner D, Peichel CL. Along the speciation continuum in sticklebacks. J Fish Biol. 2009;75:2000–36.

    Article  CAS  PubMed  Google Scholar 

  4. Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2009;24:145–56.

    Article  PubMed  Google Scholar 

  5. Peccoud J, Ollivier A, Plantegenest M, Simon J-C. A continuum of genetic divergence from sympatric host races to species in the pea aphid complex. Proc Natl Acad Sci U S A. 2009;106:7495–500.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Merrill RM, Gompert Z, Dembeck LM, Kronfrost MR, McMillan WO, Jiggins CD. Mate preference across the speciation continuum in a clade of mimetic butterflies. Evolution. 2011;65:1489–500.

    Article  PubMed  Google Scholar 

  7. Nosil P. Ecological Speciation. Oxford: Oxford University Press; 2012.

    Book  Google Scholar 

  8. Langerhans RB, Riesch R. Speciation by selection: A framework for understanding ecology’s role in speciation. Curr Zool. 2013;59:31–52.

    Article  Google Scholar 

  9. Plath M, Pfenninger M, Lerp H, Riesch R, Eschenbrenner C, Slattery PA, et al. Genetic differentiation and selection against migrants in evolutionarily replicated extreme environments. Evolution. 2013;67:2647–61.

    Article  PubMed  Google Scholar 

  10. Powell THQ, Hood GR, Murphy MO, Heilveil JS, Berlocher SH, Nosil P, et al. Genetic divergence along the speciation continuum: the transition from host races to species in Rhagoletis (Diptera: Tephritidae). Evolution. 2013;67:2561–76.

    Article  PubMed  Google Scholar 

  11. Taylor EB, Boughman JW, Groenenboom M, Sniatynski M, Schluter D, Gow JL. Speciation in reverse: morphological and genetic evidence of the collapse of a three-spined stickleback (Gasterosteus aculeatus) species pair. Mol Ecol. 2006;15:343–55.

    Article  CAS  PubMed  Google Scholar 

  12. Seehausen O, Takimoto G, Roy D, Jokela J. Speciation reversal and biodiversity dynamics with hybridization in changing environments. Mol Ecol. 2008;17:30–44.

    Article  PubMed  Google Scholar 

  13. Behm JE, Ives AR, Boughman JW. Breakdown in postmating isolation and the collapse of a species pair through hybridization. Am Nat. 2010;175:11–26.

    Article  PubMed  Google Scholar 

  14. Culumber ZW, Fisher HS, Tobler M, Mateos M, Barber PH, Sorenson MD, et al. Replicated hybrid zones of Xiphophorus swordtails along an elevational gradient. Mol Ecol. 2011;20:342–56.

    Article  CAS  PubMed  Google Scholar 

  15. Vonlanthen P, Bittner D, Hudson A, Young K, Muller R, Lundsgaard-Hansen B, et al. Eutrophication causes speciation reversal in whitefish adaptive radiations. Nature. 2012;482:357–62.

    Article  CAS  PubMed  Google Scholar 

  16. Coyne JA, Orr HA. Speciation. Sunderland: Sinauer Associates, Inc. Publishers; 2004.

    Google Scholar 

  17. Harvey PH, Pagel MD. The Comparative Method in Evolutionary Biology. Oxford: Oxford University Press; 1991.

    Google Scholar 

  18. Schluter D. The Ecology of Adaptive Radiation. Oxford: Oxford University Press; 2000. p. 2000.

    Google Scholar 

  19. Losos JB. Convergence, adaptation, and constraint. Evolution. 2011;65:1827–40.

    Article  PubMed  Google Scholar 

  20. Nosil P, Feder JL. Genome evolution and speciation: towards quantitative descriptions of pattern and process. Evolution. 2013;67:2461–7.

    Article  PubMed  Google Scholar 

  21. Reznick DN, Bryga HA. Life-history evolution in guppies (Poecilia reticulata: Poeciliidae). V. Genetic basis of parallelism in life histories. Am Nat. 1996;147:339–59.

    Article  Google Scholar 

  22. Hulsey CD, Roberts RJ, Lin ASP, Guldberg R, Streelman JT. Convergence in a mechanically complex phenotype: Detecting structural adaptations for crushing in cichlid fish. Evolution. 2008;62:1587–99.

    Article  PubMed  Google Scholar 

  23. Elmer KR, Kusche H, Lehtonen TK, Meyer A. Local variation and parallel evolution: morphological and genetic diversity across a species complex of neotropical crater lake cichlid fishes. Phil Trans R Soc B. 2010;365:1763–82.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Riesch R, Plath M, García de León FJ, Schlupp I. Convergent life-history shifts: toxic environments result in big babies in two clades of poeciliids. Naturwissenschaften. 2010;97:133–41.

    Article  CAS  PubMed  Google Scholar 

  25. Conner J, Via S. Patterns of phenotypic and genetic correlations among morphological and life-history traits in wild radish, Raphanus raphanistrum. Evolution. 1993;47:704–11.

    Article  Google Scholar 

  26. Wesner JS, Billman EJ, Meier A, Belk MC. Morphological convergence during pregnancy among predator and nonpredator populations of the livebearing fish Brachyrhaphis rhabdophora (Teleostei: Poeciliidae). Biol J Linn Soc. 2011;104:386–92.

    Article  Google Scholar 

  27. Fišer C, Zagmajster M, Zakšek V. Coevolution of life history traits and morphology in female subterranean amphipods. Oikos. 2013;122:770–8.

    Article  Google Scholar 

  28. Ghalambor CK, Walker JA, Reznick DN. Multi-trait selection, adaptation, and constraints on the evolution of burst swimming performance. Integr Comp Biol. 2003;43:431–8.

    Article  PubMed  Google Scholar 

  29. Walsh B, Blows MW. Abundant genetic variation + strong selection = multivariate genetic constraints: a geometric view of adaptation. Annu Rev Ecol Evol Syst. 2009;40:41–59.

    Article  Google Scholar 

  30. Cooper CE, Brown GC. The inhibition of mitochondrial cytochrome oxidase by the gases carbon monoxide, nitric oxide, hydrogen cyanide and hydrogen sulfide: chemical mechanism and physiological significance. J Bioenerg Biomembr. 2008;40:533–9.

    Article  CAS  PubMed  Google Scholar 

  31. Riesch R, Tobler M, Plath M. Hydrogen sulfide-toxic habitats. In: Riesch R, Tobler M, Plath M, editors. Extremophile Fishes. Heidelberg: Springer International Publishing; 2015. p. 137–59.

    Google Scholar 

  32. Tobler M, DeWitt TJ, Schlupp I, García de León FJ, Herrmann R, Feulner PGD, et al. Toxic hydrogen sulfide and dark caves: Phenotypic and genetic divergence across two abiotic environmental gradients in Poecilia mexicana. Evolution. 2008;62:2643–59.

    Article  PubMed  Google Scholar 

  33. Tobler M, Riesch R, Tobler CM, Schulz-Mirbach T, Plath M. Natural and sexual selection against immigrants maintains differentiation among micro-allopatric populations. J Evol Biol. 2009;22:2298–304.

    Article  CAS  PubMed  Google Scholar 

  34. Tobler M, Palacios M, Chapman LJ, Mitrofanov I, Bierbach D, Plath M, et al. Evolution in extreme environments: replicated phenotypic differentiation in livebearing fish inhabiting sulfidic springs. Evolution. 2011;65:2213–28.

    Article  PubMed  Google Scholar 

  35. Plath M, Hauswaldt JS, Moll K, Tobler M, de León FJ G, Schlupp I, et al. Local adaptation and pronounced genetic differentiation in an extremophile fish, Poecilia mexicana, inhabiting a Mexican cave with toxic hydrogen sulphide. Mol Ecol. 2007;16:967–76.

    Article  CAS  PubMed  Google Scholar 

  36. Plath M, Hermann B, Schröder C, Riesch R, Tobler M, García de León FJ, et al. Locally adapted fish populations maintain small-scale genetic differentiation despite perturbation by a catastrophic flood event. BMC Evol Biol. 2010;10:256.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Riesch R, Plath M, Schlupp I. Toxic hydrogen sulfide and dark caves: life-history adaptations in a livebearing fish (Poecilia mexicana, Poeciliidae). Ecology. 2010;91:1494–505.

    Article  PubMed  Google Scholar 

  38. Riesch R, Plath M, Schlupp I. Toxic hydrogen sulphide and dark caves: pronounced male life-history divergence among locally adapted Poecilia mexicana (Poeciliidae). J Evol Biol. 2011;24:596–606.

    Article  CAS  PubMed  Google Scholar 

  39. Riesch R, Schlupp I, Langerhans RB, Plath M. Shared and unique patterns of embryo development in extremophile poeciliids. PLoS One. 2011;6, e27377.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Palacios M, Arias–Rodriguez L, Plath M, Eifert C, Lerp H, Lamboj A, et al. The rediscovery of a long described species reveals additional complexity in speciation patterns of poeciliid fishes in sulfide springs. PLoS One. 2013;8, e71069.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Riesch R, Plath M, Schlupp I, Tobler M, Langerhans RB. Colonisation of toxic environments drives predictable life-history evolution in livebearing fishes (Poeciliidae). Ecol Lett. 2014;17:65–71.

    Article  PubMed  Google Scholar 

  42. Powell E. Oxygen, sulfide and diffusion: Why thiobiotic meiofauna must be sulfide-insensitive first-order respirers. J Mar Res. 1989;47:887–932.

    Article  CAS  Google Scholar 

  43. Jahn A, Janas U, Theede H, Szaniawska A. Significance of body size in sulphide detoxification in the Baltic clam Macoma balthica (Bivalvia, Tellinidae) in the Gulf of Gdansk. Mar Ecol Prog Ser. 1997;154:175–83.

    Article  CAS  Google Scholar 

  44. Bashey F. Competition as a selective mechanism for larger offspring size in guppies. Oikos. 2008;117:104–13.

    Article  Google Scholar 

  45. Passow CN, Greenway R, Arias-Rodriguez L, Jeyasingh PD, Tobler M. Reduction of energetic demands through modification of body size and routine metabolic rates in extremophile fish. Physiol Biochem Zool. 2015;88:371–83.

    Article  PubMed  Google Scholar 

  46. Smith CC, Fretwell SD. The optimal balance between size and number of offspring. Am Nat. 1974;108:499–506.

    Article  Google Scholar 

  47. Tobler M, Hastings L. Convergent patterns of body shape differentiation in four different clades of poeciliid fishes inhabiting sulfide springs. Evol Biol. 2011;38:412–21.

    Article  Google Scholar 

  48. Slattery P, Eschenbrenner C, Arias-Rodriguez L, Streit B, Bierbach D, Riesch R, et al. Twelve new microsatellite loci for the sulphur molly (Poecilia sulphuraria) and the related Atlantic molly (P. mexicana). Cons Genet Resour. 2012;4:935–7.

    Article  Google Scholar 

  49. Rosen DE, Bailey RM. The poeciliid fishes (Cyprinodontiformes), their structure, and zoogeography and systematics. Bull Am Mus Nat Hist. 1963;126:1–176.

    Google Scholar 

  50. Slatkin M. Gene flow and the geographic structure of natural populations. Science. 1987;236:787–92.

    Article  CAS  PubMed  Google Scholar 

  51. Hendry AP, Day T, Taylor EB. Population mixing and the adaptive divergence of quantitative traits in discrete populations: a theoretical framework for empirical tests. Evolution. 2001;55:459–66.

    Article  CAS  PubMed  Google Scholar 

  52. Kaeuffer R, Peichel CL, Bolnick DI, Hendry AP. Parallel and nonparallel aspects of ecological, phenotypic, and genetic divergence across replicate population pairs of lake and stream stickleback. Evolution. 2012;66:402–18.

    Article  PubMed  Google Scholar 

  53. Conner WH, Day Jr JW, Baumann RH, Randall JM. Influence of hurricanes on coastal ecosystems along the northern Gulf of Mexico. Wetlands Ecol Manage. 1989;1:45–56.

    Article  Google Scholar 

  54. Landsea CW, Pieke Jr RA, Mestas-Nuñez AM, Knaff JA. Atlantic Basin hurricanes: Indices of climatic changes. Clim Change. 1999;42:89–129.

    Article  Google Scholar 

  55. Tobler M, Scharnweber K, Greenway R, Passow CN, Arias-Rodriguez L, Garcia de León FJ. Convergent changes in the trophic ecology of extremophile fish occurring along replicated environmental gradients. Freshw Biol. 2015;60:768–80.

    Article  Google Scholar 

  56. Einum S, Fleming IA. Highly fecund mothers sacrifice offspring survival to maximize fitness. Nature. 2000;405:565–7.

    Article  CAS  PubMed  Google Scholar 

  57. Langerhans RB. Morphology, performance, fitness: functional insight into a post-Pleistocene radiation of mosquitofish. Biol Lett. 2009;5:488–91.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Langerhans RB. Predicting evolution with generalized models of divergent selection: a case study with poeciliid fish. Integr Comp Biol. 2010;50:1167–84.

    Article  PubMed  Google Scholar 

  59. Riesch R, Oranth A, Dzienko J, Karau N, Schießl A, Stadler S, et al. Extreme habitats are not refuges: poeciliids suffer from increased aerial predation risk in sulfidic, southern Mexican habitats. Biol J Linn Soc. 2010;101:417–26.

    Article  Google Scholar 

  60. Riesch R, Plath M, Schlupp I. The offspring size/fecundity trade-off and female fitness in the Atlantic molly (Poecilia mexicana, Poeciliidae). Environ Biol Fish. 2012;94:457–63.

    Article  Google Scholar 

  61. Pfenninger M, Lerp H, Tobler M, Passow C, Kelley JL, Funke E, et al. Parallel evolution of cox genes in H2S-tolerant fish as key adaptation to a toxic environment. Nat Commun. 2014;5:3873.

    Article  CAS  PubMed  Google Scholar 

  62. Tobler M, Riesch R, de León FJ G, Schlupp I, Plath M. Two endemic and endangered fishes, Poecilia sulphuraria (Alvarez, 1948) and Gambusia eurystoma Miller, 1975 (Poeciliidae, Teleostei) as only survivors in a small sulphidic habitat. J Fish Biol. 2008;72:523–33.

    Article  Google Scholar 

  63. Lydeard C, Wooten MC, Meyer A. Cytochrome b sequence variation and a molecular phylogeny of the live-bearing fish genus Gambusia (Cyprinodontiformes: Poeciliidae). Can J Zool. 1995;73:213–27.

    Article  CAS  Google Scholar 

  64. Reznick DN, Endler JA. The impact of predation on life history evolution in Trinidadian guppies (Poecilia reticulata). Evolution. 1982;36:160–77.

    Article  Google Scholar 

  65. Andersson M. Sexual selection. Princeton: Princeton University Press; 1994.

    Google Scholar 

  66. Edward DA, Chapman T. The evolution and significance of male mate choice. Trends Ecol Evol. 2011;26:647–54.

    Article  PubMed  Google Scholar 

  67. Rohlf F. tpsDig. 2010. Available from Accessed 30 June 2010.

  68. Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53:983–97.

    Article  CAS  PubMed  Google Scholar 

  69. Riesch R, Martin RA, Langerhans RB. Predation’s role in life-history evolution of a livebearing fish and a test of the Trexler-DeAngelis model of maternal provisioning. Am Nat. 2013;181:78–93.

    Article  PubMed  Google Scholar 

  70. Martin RA, Riesch R, Heinen-Kay JL, Langerhans RB. Evolution of male coloration during a post-Pleistocene radiation of Bahamas mosquitofish (Gambusia hubbsi). Evolution. 2014;68:397–411.

    Article  PubMed  Google Scholar 

  71. Langerhans RB. Trade-off between steady and unsteady swimming underlies predator-driven divergence in Gambusia affinis. J Evol Biol. 2009;22:1057–75.

    Article  CAS  PubMed  Google Scholar 

  72. Rohlf F. tpsRegr. 2009. Available from Accessed 30 June 2010.

  73. Spencer CC, Chlan CA, Neigel JE, Scribner KT, Wooten MC, Leberg PL. Polymorphic microsatellite markers in the western mosquitofish, Gambusia affinis. Mol Ecol. 1999;8:157–68.

    CAS  PubMed  Google Scholar 

  74. Purcell KM, Lance SL, Jones KL, Stockwell CA. Ten novel microsatellite markers for the western mosquitofish Gambusia affinis. Conserv Genet Resour. 2011;3:361–3.

    Article  Google Scholar 

  75. Excoffier L, Lischer H. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.

    Article  PubMed  Google Scholar 

  76. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  77. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol Ecol. 2005;14:2611–20.

    Article  CAS  PubMed  Google Scholar 

  78. Earl DA, von Holdt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.

    Article  Google Scholar 

  79. Goudet J. 2001. FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3;

  80. Riesch R, Tobler M, Lerp H, Jourdan J, Doumas T, Nosil P, Langerhans RB, Plath M. Data from: Extremophile Poeciliidae. multivariate insights into the complexity of speciation along replicated ecological gradients. Dryad Digital Repository.

Download references


We would like to thank L. Arias-Rodriguez, I. Schlupp, R. A. Martin, S. E. Diamond, and J. L. Heinen-Kay for help in the field.


Financial support came from the Erwin-Riesch Stiftung and W. M. Keck Center for Behavioral Biology (to RR), the research funding program "LOEWE − Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz" of Hesse's Ministry of Higher Education, Research, and the Arts (to MP), and the National Science Foundation (IOS-1121832, IOS-1463720, and IOS-1557860 to MT, DEB-0842364 to RBL).

Availability of data and materials

The data sets supporting the results of this article are available in the repository, [80].

Authors’ contributions

RR, RBL, PN, and MP conceived the project. RR, MT, HL, JJ, TD, RBL, and MP collected and analysed the data. All authors interpreted the data. RR, RBL and MP wrote the initial manuscript and all authors contributed to further writing and revisions. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Collection of fishes was conducted under authorization by the appropriate government agencies, which kindly provided collection permits (Mexican government: CONAPESCA: DGOPA/16986/191205/8101, DGOPA/02232/230706/1079, DGOPA.06192. 240608-1562, and SGPA/DGVS/04751/08; the Municipal of Tacotalpa: SM/1133/208; the United States National Park Service Chickasaw NRA: CHIC-2007-SCI-0001), and the study is compliant with all relevant laws for the ethical treatment of animals in scientific research.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Rüdiger Riesch.

Additional file

Additional file 1:

OSM 1: Additional information on Gambusia collection sites. (DOCX 2179 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Riesch, R., Tobler, M., Lerp, H. et al. Extremophile Poeciliidae: multivariate insights into the complexity of speciation along replicated ecological gradients. BMC Evol Biol 16, 136 (2016).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: