Figures
Abstract
The process of ecological speciation drives the evolution of locally adapted and reproductively isolated populations in response to divergent natural selection. In Southern Mexico, several lineages of the freshwater fish species of the genus Poecilia have independently colonized toxic, hydrogen sulfide-rich springs. Even though ecological speciation processes are increasingly well understood in this system, aligning the taxonomy of these fish with evolutionary processes has lagged behind. While some sulfide spring populations are classified as ecotypes of Poecilia mexicana, others, like P. sulphuraria, have been described as highly endemic species. Our study particularly focused on elucidating the taxonomy of the long described sulfide spring endemic, Poecilia thermalis Steindachner 1863, and investigates if similar evolutionary patterns of phenotypic trait divergence and reproductive isolation are present as observed in other sulfidic species of Poecilia. We applied a geometric morphometric approach to assess body shape similarity to other sulfidic and non-sulfidic fish of the genus Poecilia. We also conducted phylogenetic and population genetic analyses to establish the phylogenetic relationships of P. thermalis and used a population genetic approach to determine levels of gene flow among Poecilia from sulfidic and non-sulfidic sites. Our results indicate that P. thermalis' body shape has evolved in convergence with other sulfide spring populations in the genus. Phylogenetic analyses placed P. thermalis as most closely related to one population of P. sulphuraria, and population genetic analyses demonstrated that P. thermalis is genetically isolated from both P. mexicana ecotypes and P. sulphuraria. Based on these findings, we make taxonomic recommendations for P. thermalis. Overall, our study verifies the role of hydrogen sulfide as a main factor shaping convergent, phenotypic evolution and the emergence of reproductive isolation between Poecilia populations residing in adjacent sulfidic and non-sulfidic environments.
Citation: Palacios M, Arias-Rodriguez L, Plath M, Eifert C, Lerp H, Lamboj A, et al. (2013) The Rediscovery of a Long Described Species Reveals Additional Complexity in Speciation Patterns of Poeciliid Fishes in Sulfide Springs. PLoS ONE 8(8): e71069. https://doi.org/10.1371/journal.pone.0071069
Editor: Norman Johnson, University of Massachusetts, United States of America
Received: April 24, 2013; Accepted: July 1, 2013; Published: August 16, 2013
Copyright: © 2013 Palacios et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by NSF-IOS (Grant number: IOS-1121832, http://www.nsf.gov/div/index.jsp?div=ios); and LOEWE – Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz” of Hesse's Ministry of Higher Education, Research, and the Arts (Grant number: Bik-F Projektbereich C2, http://verwaltung.hessen.de/irj/HMWK_Internet?uid=fa560c0b-ed11-9311-1010-43bf5aa60dfa). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Divergent natural selection, often mediated by environmental variation, is a key driver of phenotypic evolution. Its effects can lead to the emergence of locally adapted populations that exhibit unique traits and occupy habitats with distinct combinations of environmental factors [1]. Depending on the strength of selection and rates of gene flow, such local adaptation may also cause the emergence of reproductive isolating barriers among diverging populations as a byproduct, a process known as ecological speciation [2]–[4]. Such speciation driven by adaptation has been documented in a variety of organisms and in response to a diversity of selective forces (e.g., resource exploitation [5], [6], habitat use [7], [8], and predation [9], [10]). Current research efforts focus on elucidating the range of conditions under which local adaptation is likely to translate into reproductive isolation, particularly because speciation does not appear to be an inevitable consequence of divergent selection even when phenotypic differentiation is pronounced [11]–[13]. The fact that phenotypic divergence is not necessarily tied to the emergence of reproductive isolation, as well as the dynamic nature of species boundaries, pose a challenge for efforts in biological systematics and taxonomy, which attempt to describe and organize biological diversity [14].
Historically, taxonomists catalogued species solely based on morphological trait variation (i.e., they applied a morphological species concept) without considering the underlying mechanisms contributing to phenotypic variation and speciation. Darwin [15] himself conceded that the recognition of species was often left to the opinion, experience, and expertise of naturalists. Confounding evolutionary processes and taxonomy has led to much confusion for both taxonomists and evolutionary biologists [16]. While the Modern Synthesis consolidated the definition of a biological species to being based on reproductive isolation [17], [18], a wide variety of other species concepts is currently being applied and newly developed depending on the objective at hand [19]. Consequently, there is a perpetual struggle to align the works of earlier biologists with a more modern understanding of phenotypic evolution and speciation [19], particularly in diverse and taxonomically difficult groups with a wealth of available species names, many of which currently are considered junior synonyms of older epithets [20]. Reexamination of long described taxa is critical to accurately assess biodiversity and align taxonomy with evolutionary processes leading to diversity, because phenotypic variation due to developmental plasticity [21], genetically based intraspecific polymorphisms [22], and large-scale geographic trait variation are not necessarily tied to reproductive isolation.
The genus Poecilia, which is part of the livebearer family Poeciliidae, represents an excellent example of the difficulties in aligning evolutionary processes and taxonomy. Poecilia is a diverse group of freshwater fish species that are distributed from the southeastern United States and Middle America to parts of South America and the Greater Antilles [23]. Taxonomic confusion is particularly prevalent in the P. sphenops (short fin molly) species group, which occurs from northern Mexico to Venezuela [24]. On one hand, these fish show tremendous variability in phenotypic traits, which has lead to the description of numerous species (many of which are currently considered synonyms) [25]. Despite clear phylogenetic structuring [26] and a long list of available names [25], species designation often remains unclear. This is predominantly caused by pronounced eco-morphological and geographic variation [27], intra-specific trait variability that appears to regularly exceed the inter-specific differences [23], [28], [29], uncertainty about whether and how morphological differences are tied to reproductive isolation [30], [31], potential hybridization and introgression among lineages [26], [32], [33], and sometimes unclear type localities of available names [25]. On the other hand, there is some well documented cases of phenotypically divergent and reproductively isolated species; yet these lineages often do not appear to be phylogenetically distinct, presumably because divergence occurred relatively recently [34]. Therefore, evaluation of species in this group in an evolutionary context is crucial for the resolution of taxonomy and for a better understanding of the historical and current processes shaping the species' evolution.
Here, we attempt to clarify the status of one such long described species, Poecilia thermalis Steindachner 1863 (Fig. 1), which has had a vivid taxonomic history (Table 1). The species was originally described based on specimens that C. B. Heller collected in a sulfidic spring (La Esperanza) located in the Ixtapangajoya river drainage of Chiapas, Mexico, in 1848 [35]. The species description particularly emphasized the large head size in the available specimens. Shortly after the species description, Günther [36] thought to recognize P. thermalis in samples from El Salvador, but this has been largely viewed as a misidentification [37]. Subsequently, the species was subjected to various nomenclatural re-assignments by taxonomists (i.e., it has predominantly been viewed as a synonym to P. sphenops; see Table 1). Interestingly, all prior nomenclatural acts revolving around this species have either been based solely on the consultation of the misidentified specimens [23], [38], [39] or on information given in the species description (without actually examining specimens [24], [37]). Because the species has not been collected since 1848, we revisited the type locality of P. thermalis based on descriptions in Heller's [40] autobiographical account of his travels in southern Mexico to investigate the status of this enigmatic species. Our interest was fueled by the sulfidic nature of P. thermalis' habitat, which was emphasized both in the species description and Heller's field accounts [35], [40]. In drainages adjacent to the Río Ixtapangajoya (namely the Tacotalpa and Puyacatengo drainages to the east, and the Pichucalco drainage to the west), evolutionarily independent Poecilia lineages have colonized springs with high concentrations of toxic hydrogen sulfide (H2S) [34]. These sulfide spring inhabitants are phenotypically distinct from the closely related P. mexicana in non-sulfidic environments within the same drainage and are characterized by morphological, physiological, behavioral, and life history adaptations that show strong signals of convergent evolution across drainages [34]. In particular, sulfide spring fishes are characterized by enlarged heads and correlated increases in gill surface area, which facilitates oxygen acquisition in hypoxic sulfide spring environments and directly affects survival [34], [41], as well as physiological and biochemical adaptations that reduce the impacts of sulfide toxicity [42]. In conjunction with adaptive trait divergence, the sulfide spring populations are reproductively isolated from adjacent populations from non-sulfidic waters despite small geographic distances (in some instances <100 meters) and a lack of physical barriers that would prevent fish migration [43], [44]. Reproductive isolation appears to be mediated particularly by natural and sexual selection against immigrants [43]. The taxonomic status of sulfide spring populations varies across drainages; sulfide spring residents in the Tacotalpa and Puyacatengo drainages are considered ecotypes of the widespread P. mexicana despite clear morphological differences and strong reproductive isolation, while sulfide spring residents in the Pichucalco drainage have been described as a distinct species, Poecilia sulphuraria [45].
A.-F. Representative specimens of Poecilia thermalis (males on top, females below). A. & B. represent artistic depictions of the species from the original description by Steindachner [35]. C.-F. are photos of freshly collected specimens of the P. thermalis from the large La Esperanza spring in 2012. G. The type locality of P. thermalis at La Esperanza (large spring), Chiapas, Mexico, a typical sulfidic spring habitat. H. A second, smaller sulfide spring (La Esperanza small spring) in close proximity to the type locality, which is also inhabited by P. thermalis.
Consequently, this study investigated the re-discovered sulfide spring species Poecilia thermalis to examine whether it shows similar evolutionary patterns (i.e., phenotypic trait divergence and reproductive isolation) as other sulfidic populations in the region and to shed light on its taxonomy. Specifically, we used morphological, phylogenetic, and population genetic approaches to address three major questions: (1) How do specimens from the type locality of Poecilia thermalis phenotypically compare to other Poecilia populations from sulfidic and non-sulfidic spring habitats in the region? Using a geometric morphometric approach, we tested for potential morphological convergence between P. thermalis and other sulfide spring fish in southern Mexico. We also explored the similarity of body shape between historical samples of P. thermalis [35] to current populations of sulfidic and non-sulfidic spring fish, including recently collected P. thermalis from the type locality. (2) What is the phylogenetic relationship of P. thermalis to other mollies? Based on mitochondrial and nuclear markers from a broad taxonomic sampling of Poecilia, we elucidated the phylogenetic position of P. thermalis. We were particularly interested in determining whether the species represents a unique sulfide-adapted lineage within Poecilia. (3) Is P. thermalis genetically isolated from adjacent Poecilia populations? Sulfide spring populations of Poecilia consistently exhibit a high degree of reproductive isolation from non-sulfidic populations despite the small spatial distance and a lack of migratory barriers [46]–48. Using a population genetic approach based on microsatellites, we quantified gene flow between the endemic sulfide-adapted species of P. thermalis and P. sulphuraria, and P. mexicana populations from adjacent non-sulfidic habitats.
Materials and Methods
Ethics Statement
All procedures conducted for this study were approved by the Institutional Animal Care and Use Committee at Oklahoma State University (ACUP: AS10-15) and permits issued by the Municipio de Tacotalpa-Tabasco (DFET/23/2011), as well as the Mexican Federal Agencies SEMARNAT (SGPA/DGVS/04315/11 for Poecilia sulphuraria) and CONAPESCA (DGOPA.09004.041111.3088 for Poecilia sp.). Details about the selection of focal populations are given below for each section separately.
Study sites and sample collection
Our study area lies in the foothills of the Sierra Madre de Chiapas in the northeastern part of the state of Chiapas, where the mountains meet the wide floodplains of the state of Tabasco. Here, four tributaries of the Río Grijalva, the Tacotalpa, Puyacatengo, Ixtapangajoya, and Pichucalco (from east to west), provide a system of naturally replicated non-sulfidic and adjacent sulfidic habitats (see Figure 2). The presence of high concentrations of hydrogen sulfide (H2S) in the sulfidic springs is a result of a nearby, active volcano, El Chinchón, and bacterial activity [49]–[51]. Besides the presence of H2S, the sulfidic spring environments differ from non-sulfidic habitats in a variety of environmental factors, including reduced oxygen concentration, structural habitat differences, reduced species richness, and reduced photoautotrophic primary production [34], [52], [53]. In all cases, sulfidic habitats drain directly into non-sulfidic habitats in the same drainage, and there are no major physical barriers preventing movement of fish, especially during high flow periods. The distances between sulfide spring habitats to the closest freshwater habitat range from <50 m to about 500 m.
Sampling of P. thermalis, Poecilia mexicana mexicana, and P. sulphuraria in southern Mexico. The colors represent sulfidic (yellow) and non-sulfidic (blue) sites and the numbers represent localities as described in Table 2 of the main manuscript. For orientation purposes, we included black lines representing major roads and gray areas representing major towns in the region.
While abiotic environmental parameters of other sites investigated here have been published in detail before [34], the type locality of P. thermalis at ‘La Esperanza’ has not been resampled since Heller's expedition. Exploration of the general area in May and June of 2012 revealed not one, but two proximate sulfide springs, a larger one (site 7) and a smaller one (site 8) separated by a freshwater tributary (site 9) of the Río Ixtapangajoya (Fig. 2). Based on descriptions in his autobiography, it is evident that Heller collected the type specimens from the larger spring (hence this larger spring should be considered as the type locality), which represents a cluster of relatively high discharge sulfide springs that forms a distinct tributary flowing over gravelly bottom for about 200 meters before merging with the Río Ixtapangajoya over a small drop. In this habitat, we measured the following water parameters (N = 4 for all measurements): sulfide concentration: 216±38 μM; dissolved oxygen: 1.31±0.39 mg/l; pH: 6.9±0.0; specific conductivity: 4.132±0.083 mS/cm; temperature: 28.7±0.9°C. The smaller spring represents a group of lower discharge springs in a swampy area with dense reeds. The site is merely a shallow pool (about 15×20 m) with a narrow outflow that – at least during our visit at the end of the dry season – eventually disappeared into pasture grounds (i.e., there was no direct, permanent connection to the Río Ixtapangajoya less than 150 m away). Here, we measured the following water parameters (N = 4 for all measurements): sulfide concentration: 41±3 μM; dissolved oxygen: 1.67±0.21 mg/l; pH, 7.0±0.1; specific conductivity: 2.979±0.071 mS/cm; temperature: 27.2±0.3°C. Overall, the physiochemical conditions in the two La Esperanza springs aligned well with data collected over multiple years in sulfide springs of other drainages [34].
All specimens for this study were collected using seines, euthanized with buffered MS222 immediately after capture, and fixed in a 10% formaldehyde solution for geometric morphometric analyses. All specimens are housed in the Department of Zoology, Oklahoma State University. In addition, we took fin clips (right pectoral fin) that were preserved in 95% ethanol and stored at 4°C for molecular analyses.
Geometric morphometrics
The specimens sampled for the geometric morphometric analyses broadly covered the study region, including all currently known sulfide springs inhabited by Poecilia species. Our sampling included the southern subspecies of P. mexicana mexicana from a variety of non-sulfidic habitats in all drainages, sulfidic ecotypes of P. m. mexicana in the Tacotalpa and Puyacatengo drainages as well as the previously described sulfide spring endemics P. thermalis (Ixtapangajoya drainage) and P. sulphuraria (Pichucalco drainage; see Table 2 for an overview). We analyzed a total of 1099 specimens from the four drainages for a total of 23 localities.
We conducted a geometric morphometric analysis of body shape both from the lateral and dorsal view, as body shape has been shown to be a reliable indicator for convergent evolution in sulfide spring environments [34], [54]. For all specimens, lateral and dorsal photographs were taken using a Nikon D90 camera mounted on a copy stand. We digitized 16 lateral and 9 dorsal landmark points using tpsDig2 [55] (see Figure S1A, 1B in File S1 for details on landmark locations). We analyzed lateral and dorsal landmarks separately and performed a geometric morphometric analysis based on the coordinates of the digitized landmarks [56]. Landmark coordinates were aligned using least-square superimposition as implemented in the program tpsRelw [57] to remove effects of translation, rotation, and scale. Based on the aligned coordinates, we calculated centroid size and partial warp scores with uniform components (weight matrix) for each individual. Unless otherwise stated, all statistical analyses were performed using SPSS 20 (IBM Inc.). All raw data used for the analyses described below are archived on http://datadryad.org under the DOI number associated with this publication.
The weight matrices obtained from the geometric morphometric analyses were first subjected to principal components analyses (PCA) using a covariance matrix to reduce data dimensionality. We retained 9 PC axes with an eigenvalue greater than 1 for the lateral dataset (explaining >95% of variation) and 8 PC axes for the dorsal dataset (explaining >96% of variation). Individual PC axis scores were used as dependent variables in multivariate analyses of covariance (MANCOVA). Assumptions of multivariate normal error and homogeneity of variances and co-variances were met for all analyses performed. F-values were approximated using Wilks' lambda and effect strengths by use of partial eta squared (ηp2). We also calculated the relative variance as the partial variance for a given term divided by the maximum partial variance value in a model [58]. We included sex, habitat type (H2S present or not), drainage, and site (nested within the H2S × drainage interaction) as well as all interaction terms as independent variables. Centroid size was included in the models as a covariate to control for multivariate allometry.
Since random nested factors are not applicable for MANCOVAs, and the use of fixed effects can inflate type I error rates when nested terms are significant, we also analyzed shape variation using a mixed-model nested analysis of covariance (ANCOVA) [59]. To do so, we calculated divergence scores for each individual along the sulfide/non-sulfide gradient based on a divergence vector as defined by Langerhans [59]. Individual divergence scores were used as dependent variables in ANCOVAs with the same model structure as outlined above, except that site was designated a random factor to account for the fact that only a random subset of sites where Poecilia occurs was analyzed for this study. Shape variation along the first two PC axes and along the sulfide/non-sulfide divergence axes was visualized with thin-plate spline transformation grids using tpsRegr [60].
To fully examine the multidimensional affinities of different Poecilia populations relative to each other, including information from both the lateral and the dorsal projection, weight matrices for both projections were combined and subjected to a principal components analysis, from which we retained 14 axes with an Eigenvalue >1. Population-specific estimated marginal means for each axis were calculated using a MANCOVA model as detailed above and used to create a dissimilarity matrix that was subjected to a hierarchical cluster analysis using the neighbor-joining algorithm [61].
Finally, we tested whether the type specimens of P. thermalis collected in 1848 by Heller clustered with specimens we collected from the two La Esperanza springs in 2012. To do so, lateral photographs were taken from the available syntypes in the collection of the Natural History Museum in Vienna (N = 18), and we digitized the same lateral landmarks as for all other specimens. We compared the museum specimens to the samples obtained from the two La Esperanza sulfide springs, the two most proximate non-sulfidic locations in the same drainage (tributary to Río Ixtapangajoya and Río Ixtapangajoya proper), and – given the phylogenetic affinity of P. thermalis to P. sulphuraria (see below) – to the specimens obtained from the Baños del Azufre and La Gloria sulfide springs. Landmark coordinates from said collections were aligned separately. The weight matrix was then subjected to PCA, and the effects of sex and allometry removed from the dataset by using the residuals of a preparatory MANCOVA, in which the principal component scores were used as dependent variables, centroid size as a covariate, and sex as an independent variable. We then conducted a discriminant function analysis (DFA) to elucidate whether museum specimens were classified to the Esperanza sulfide springs based on body shape data. We used a cross-validation technique where discriminant functions were generated based on the data of contemporary samples (training data), and classification probabilities of museum specimens (testing data) to any of the six populations were calculated based on the established functions.
Phylogenetic analyses
To establish the phylogenetic relationships of P. thermalis, we sequenced a set of genes in specimens from select non-sulfidic and sulfidic habitats included in the morphometric analyses (Table S1 in File S1). In addition, we broadened our taxon sampling by adding several other species of the subgenus Mollienesia (including the endemic sulfide spring species P. sulphuraria, the Southern and Northern Mexican subspecies P. m. mexicana and P. m. limantouri, as well as P. butleri, P. sphenops, P. latipinna, and P. caucana) and more distant groups in the genus Poecilia (sensu lato; including, Limia vittata, L. dominicensis, L. melanogaster, Acanthophacelus reticulata, A. wingei, Micropoecilia bifurca, M. parae, Pamphorichthys hollandi, and P. minor). The distantly related species Cnesterodon decemmaculatus and C. hypselurus were used as outgroups to root phylogenetic trees. A complete list of all taxa examined, along with locality information and GenBank Accession numbers, is provided in Supplementary Table S1 in File S1.
The total genomic DNA was extracted from ethanol-preserved fin clips with the DNeasy Blood & Tissue Kit (Qiagen, Inc., Valencia, CA) following the manufacturer's protocol. The samples were amplified for several presumably neutral genes commonly used for phylogenetic reconstruction in fishes. Focal genes included the mitochondrial cytochrome b gene (cyt b, 1,140 base pairs) with LA and HA primers [62], the mitochondrial gene NADH subunit 2 (ND2, 1,047 bp) with ND2B-L [63] and ASN [64] primers. The nuclear genes amplified included exon 3 of recombination activating gene-1 (Rag1, 1,561 bp), a portion of the 7 trans-membrane receptor region of Rhodopsin (Rh, 822 bp), and exon 1 of myosin heavy polypeptide 6 (myh6, 767 bp) with the primers and protocol following previously published PCR protocols [65], [66]. PCR products were purified with Exosap-IT enzyme reaction (GE Healthcare Bio-Sciences Corp., Piscataway, NT), directly sequenced with a dye-labeled terminator kit (Big Dye Terminator version 3.1, Applied Biosystems, Foster City, CA), and run on an ABI automated sequencer (Applied Biosystems, Foster City, CA). Sequence electrophenograms were edited with Sequencher version 4.8 (Gene Codes) and aligned with MAFFT v. 6.0 [67].
We tested for incongruence between mitochondrial (mtDNA) and nuclear (nDNA) markers to determine evidence of introgression. Given the agreement between both datasets (data not shown), we used a concatenated dataset for further analyses. We used MrModeltest version 2.3 [68] to determine the most likely model of DNA substitution among 24 candidate models on a fixed BioNJ-JC tree based on the Akaike information criterion (AIC) (Table S2 in File S1). We also compared likelihood scores between Bayes runs of 12,000,000 generations of mtDNA as a single unit, partitioned by gene and by position to determine the most informative partition. The best likelihood score was observed in the codon partition dataset and used for analyses where possible.
For maximum likelihood (ML) analyses, we used RAxML GUI version 1.0 [69], [70] run to conduct 500 Rapid Bootstrap searches followed by an ML search. We ran the complex general time reversible (GTR) + Γ (Gamma distribution for rate variation among sites) model because RAxML does not implement simpler models. We also used GARLI version 2.0 [71] to perform ML bootstrap searches (500 replicates) on the concatenated datasets under the corresponding best model selected (Table S2 in File S1). The bootstrap trees were summarized with a Sumtrees script with a 50% percent majority rule consensus parameter in DendroPy 3.10.1 [72].
Bayesian analyses were run twice independently in MrBayes version 3.2.1 [73], [74], under models of nucleotide substitution uniquely defined for the partition of each data set (Table S2 in File S1) implementing two runs with four chains under default parameters. Appropriate “burn-in” (i.e., samples discarded prior to reaching a stationary posterior distribution) was determined based on small and stable average standard deviation of the split frequencies, potential scale reduction factor close to 1 (see MrBayes manual), and stable posterior probability values examined in Tracer version 1.5 [75]. Pairwise genetic distances based on the concatenated dataset were calculated under the Kimura-2 parameter in MEGA version 5 [76] with pairwise deletion for missing data.
Population genetics
Besides the two populations of P. thermalis, samples for population genetic analyses included two proximate non-sulfidic sites within the Ixtapangajoya drainage to test for gene flow between adjacent populations from sulfidic and non-sulfidic waters. Given the phylogenetic clustering of P. thermalis with P. sulphuraria from the Pichucalco drainage (see below), we also included both known P. sulphuraria populations as well as adjacent P. mexicana samples from that drainage into our analyses (Table 2). We used 17 previously developed microsatellite markers [77], [78] to genotype a total of 272 samples and arranged the microsatellites into three multiplex reactions [43]. Data from 80 specimens (Baños del Azufre and Puente El Azufre II) were re-analyzed from a previous study [43]. All raw data used for the population genetic analyses are archived on http://datadryad.org under the DOI number associated with this publication.
We extracted DNA using the NucleoSpin®Tissue kit (Macherey-Nagel). Microsatellites were amplified with the Type-it Microsatellite PCR kit from Qiagen (Hilden, Germany). The PCR protocol included an initial denaturation step for 5∶00 min at 95°C, 30 cycles of 1∶30 min at 60°C, and 0∶30 min at 72°C, followed by a final extension step for 30∶00 min at 60°C. The 5 µl reaction mix included 2.5 µl Type-it master mix, 0.4 µl primer mix, 0.4 µl Q-solution, 0.9 µl RNase-free water, and 0.8 µl template DNA. PCR products were analyzed on a CEQ2000 sequencer (Beckman) Coulter; denaturation at 90°C for 2 min, injection at 2.0 kV for 30 s, separation at 6.0 kV for 45 min) along with the manufacturer's internal size standard. Samples were screened using Genome Lab GeTX 10.2 software (Beckman Coulter) and alleles were called manually.
We employed the software STRUCTURE 2.3.4 [79] to identify the number of genetically distinct clusters (K) and then used the method of Evanno et al. [80] and the web-based software STRUCTURE HARVESTER 0.6.93 [81] to detect the uppermost level of population differentiation. In addition, we calculated pairwise FST-values between all population pairs and conducted a Principal Component Analysis (PCA) to further examine genetic distinctiveness between populations using GenAlEx 6.5 [82].
Results
Phenotypic variation
For the lateral projection (N = 1099 individuals), body shape varied in the position of the anal fin along PC axis 1 and the head size along PC axis 2 (Fig. 3). ‘Sex’ and the ‘presence of H2S’ explained the majority of body shape variation in our dataset (Table 3). ‘Sex’ particularly accounted for variation along the first PC axis (males have a more anterior anal fin position, as this fin is modified into a copulatory organ, the gonopodium). ‘Presence of H2S’ explained variation along the second PC axis (with sulfide spring fish having larger heads than fish from non-sulfidic habitats). All other factors and the interaction terms also had significant effects on body shape, but only ‘centroid size’, ‘drainage’, ‘site’, and the interaction of ‘H2S × drainage’ explained an appreciable amount of variation (relative variance >0.1; Table 3).
Depicted are mean principal component scores along the first two principal component axes for each site for P. thermalis (yellow circles), P. sulphuraria (yellow triangles), as well as sulfidic and non-sulfidic populations of P. mexicana across the 23 study sites in southern Mexico. The thin-plate spline transformation grids represent shape variation along each principal component axis.
These general patterns, and particularly the strong differentiation between ecotypes from sulfidic and non-sulfidic springs, were confirmed in the analysis of divergence vector scores with ‘site’ being treated as a random factor (Table 4). Visualization of shape variation along the sulfide-non-sulfide divergence vector corroborated head size as the primary difference between ecotypes (Fig. 4). Inspection of divergence scores indicated that P. thermalis, just like Poecilia from the Tacotalpa and the Pichucalco drainages, exhibits a body shape typical for sulfide spring populations. In the Puyacatengo drainage, differentiation between sulfidic and non-sulfidic ecotypes was less clear-cut, with some populations exhibiting more intermediate body shapes (Fig. 4).
Depicted are the mean divergence scores (± SEM; derived from the H2S term in the MANCOVA) for each site for the three formal species (P. thermalis, P. sulphuraria, and P. mexicana) across the 23 sites in southern Mexico from sulfidic (yellow) and non-sulfidic (blue) populations including. The numbers correspond to sites as described in Table 2.
For the dorsal projection (N = 1093), analyses revealed that body shape particularly varied in head length and width, mouth width, and body width at the insertion of the pelvic fins (see Fig. 5). ‘Presence of H2S’ explained most variation in body shape, with ecotypes from sulfidic and non-sulfidic habitats particularly segregating along the first PC axis. As for the lateral projection, all other factors and the interaction terms were also significant predictors of body shape, but only ‘sex’, ‘drainage’, ‘site’, and the interaction of ‘H2S × drainage’ had relative variance >0.1 (Table 5).
Depicted are mean principal component scores along the first two principal component axes for each site for P. thermalis (yellow circles), P. sulphuraria (yellow triangles), as well as sulfidic and non-sulfidic populations of P. mexicana across the 23 study sites in southern Mexico. The thin-plate spline transformation grids represent shape variation along each principal component axis.
The strong differentiation between ecotypes in dorsal body shape was confirmed in the analysis of divergence vector scores with site being treated as a random factor (Table 6). Visualization of shape variation along the sulfide-non-sulfide divergence vector indicated that sulfide spring fish had longer heads, wider mouths, but narrower bodies (Fig. 6). Differentiation between ecotypes from sulfidic and non-sulfidic springs was highly significant for all sites and drainages investigated, including both P. thermalis populations. The only exception was site 14 (Río Puyacatengo road crossing), which exhibited an intermediate morphology (Fig. 6).
Depicted are mean divergence scores (± SEM; derived from the H2S term in the MANCOVA) for each site for the three formal species (P. thermalis, P. sulphuraria, and P. mexicana) across the 23 sites in southern Mexico from sulfidic (yellow) and non-sulfidic (blue) populations. Numbers correspond to sites as described in Table 2.
The cluster analysis based on the combined lateral and dorsal datasets grouped 7 of 9 populations from sulfidic habitats together in a discrete cluster, highlighting a strong convergent pattern of body shape evolution in sulfide springs (Fig. 7). The two notable exceptions were fish from the La Lluvia big spring and the Puyacatengo springs (both in the Puyacatengo drainage), which were nested within non-sulfidic populations (Fig. 7). Most importantly, however, the P. thermalis samples from both La Esperanza springs in the Ixtapangajoya drainage formed a cluster with the P. sulphuraria samples from the Baños del Azufre and the La Gloria springs in the Pichucalco drainage.
Colors denote sulfide-adapted (yellow) and non-adapted (blue) populations of three species, Poecilia thermalis (yellow circles), P. sulphuraria (yellow triangles), and P. mexicana. The shapes represent the drainages (diamonds- Tacotalpa, squares- Puyacatengo, circles- Ixtapangajoya, and triangles- Pichucalco) and the numbers correspond to sites as described in Table 2.
Finally, the DFA (subset of N = 397) clearly grouped the 18 type specimens, for which we were able to obtain lateral body shape data, with contemporary sulfide spring samples, not with proximate non-sulfidic samples. Nonetheless, only 72.2% of samples were assigned to the large La Esperanza spring (Fig. S2 in File S1; Table S3 in File S1), with the remaining individuals grouped either with P. sulphuraria samples from the La Gloria springs (16.7%) and Baños del Azufre (11.1%). It should be noted at this point that the sample size for historical specimens was relatively low, such that within population variation in body shape may be underestimated. The examination of additional syntypes, which were not available at the time of our study, consequently could lead to a lower overall classification success.
Phylogenetic analyses
The evolutionary relationships in our phylogenetic analyses corroborate previously observed relationships in Poecilia among the main lineages of Acanthophacelus, Micropoecilia, Limia, Pamphorichthys, and Mollienesia [65], [66]. The results also show similar patterns previously observed in the relationships among species within the subgenus Mollienesia, with the sailfins, P. latipinna and P. latipunctata, forming a monophyletic group (Fig. 8; 100% BSS; 100% BPP). They are closely related to the monophyletic the shortfin group (100% BSS; 100% BSP [34]) with an average genetic divergence of 7.6% (Table S4 in File S1). Within the shortfin mollies, there is a separation among the P. sphenops clade (P. catemaco and P. sphenops; 100% BSS; 100% BSP) and the P. mexicana clade (P. butleri, P. sulphuraria, P. thermalis, P m. mexicana, P. m. limantouri; 100% BSS; 100% BSP) with an average genetic divergence of 6.7%. Phylogenetic analyses strongly (100% BSS and 100% BPP) indicate that P. mexicana, P. sulphuraria, and P. thermalis represent a monophyletic group. However, we did not find P. sulphuraria to be monophyletic, as P. thermalis is most closely related to the P. sulphuraria from the Baños de Azufre population, a relationship that is highly supported (87% BSS; 100% BSP). Genetic divergence between P. thermalis and this population of P. sulphuraria was generally low (0.1%) and only slightly less than the divergence to the La Gloria population (0.2%; Table S4 in File S1).
Phylogenetic analyses of two mitochondrial and three nuclear genes (5337 base pairs) yielded nodal support values (in percent) represent (from top to bottom) Bayesian Posterior Probabilities, as well as RAxML, and GARLI bootstrap support values. Asterisks denote nodal support of ≥95% for all three methods. Nodes with no values present either had low values or were of little interest for this study.
Population genetic analyses
Our Bayesian clustering analysis uncovered K = 2 as the uppermost hierarchical level of population structure. Poecilia sulphuraria (sensu lato; i.e, including both the Baños del Azufre and the La Gloria population) and P. thermalis together formed one genetic cluster that was distinct from all P. mexicana populations (Fig. 9). The only exception was one animal caught in the small Esperanza spring, which was assigned to P. mexicana, not P. thermalis. Another peak for ln P (X|K) – i.e., the second most likely level of population structure according to Evanno et al. [80] – was obtained for K = 7. In addition to detecting population genetic structure within and between drainages in P. mexicana a clear separation between P. thermalis and P. sulphuraria (s.l.) became apparent (Fig. 9), indicating low recurrent gene flow between them. Pairwise FST-values (Table S5 in File S1) revealed significant genetic differentiation between populations of P. thermalis and P. sulphuraria (0.086–0.103) as well as between populations within each species (P. thermalis, 0.031; P. sulphuraria, 0.082). Further support for genetic differentiation between the two species was obtained from the individual-based PCA (Fig. S3 in File S1).
The top panel is a bar plot showing the assignment scores of individuals by STRUCTURE with K = 2 with yellow representing P. thermalis and P. sulphuraria as a cluster from the Puyacatengo and Pichucalco sulfidic drainages and blue representing P. m. mexicana from non-sulfidic sites in the same drainage. The bottom panel is a bar plot showing the assignment scores clustering at K = 7, the second most likely number of distinct groups.
Discussion
The re-discovered species, Poecilia thermalis, is a highly endemic sulfide spring fish found to inhabit two proximate sulfidic springs in the Ixtapagajoya drainage of southern Mexico. Morphological features exhibited by P. thermalis closely resemble those of other sulfide spring fishes, independent of their phylogenetic relationship, highlighting strong patterns of convergent morphological evolution in sulfidic environments. Phylogenetic analyses placed P. thermalis as sister taxon to one (Baños del Azufre) of the two populations of the sulfidic spring endemic P. sulphuraria [45] from the Pichucalco drainage with low – albeit significant – genetic distance. In addition, population genetic analyses detected little or no current gene flow between P. thermalis and P. m. mexicana from adjacent non-sulfidic habitats or P. sulphuraria from sulfidic springs in the Pichucalco drainage, respectively, indicating an independent evolutionary trajectory of P. thermalis.
Morphological variation
In general, body shape analyses of sulfidic and non-sulfidic fishes in southern Mexico demonstrated a clear differentiation between habitat types across drainages, except for the Puyacatengo drainage. In the latter, phenotypic differentiation – particularly in the lateral projection – was less pronounced and more gradual, which could be driven by the higher spatial heterogeneity in the presence of H2S (see [43] for a discussion). Our results support previous analyses suggesting that fish in the Puyacatengo drainage were not as differentiated along a speciation continuum as sulfide spring fish from other drainages, because either they colonized sulfide springs more recently or gene flow across habitat types constrains phenotypic divergence [43]. Despite the strong patterns of convergence, we also found significant differences in body shape between sulfidic fish among drainages, indicating that there is both convergent and non-convergent aspects trait differentiation in response to sulfide exposure.
The phenotypic analyses of lateral body shape indicated that fishes from sulfidic and non-sulfidic habitats primarily differ in head size, with populations from sulfidic habitats having significantly larger heads. As such, our study validated previous findings from the genus Poecilia [34] and other poeciliids [54] with the most comprehensive sampling of sulfide spring mollies to date. Previous studies have shown that head size is positively correlated gill surface area, which is adaptive in sulfidic environments because sulfide springs exhibit low oxygen concentrations and sulfide detoxification requires additional oxygen [34]. Modification of respiratory morphological traits in conjunction with changes in respiratory behavior represent a critical adaptation mediating survival in the sulfidic and hypoxic environments [41], [83].
Our study is the first to examine body shape variation in the dorsal projection and found sulfide spring fishes to exhibit longer and wider heads, wider mouths, and narrower bodies. Sulfidic spring fishes are known to rely on aquatic surface respiration (ASR), i.e., they skim the water from the air-water interface (with higher dissolved oxygen concentrations) using their gills [41], [83]. Wider heads and mouths likely are adaptive, because they maximize the uptake of surface water as reported in neotropical characids [84] and serrasalmids [85], [86], which exhibit temporary dermal swellings of the lower jaw allowing for an increased efficiency of ASR when exposed to hypoxic conditions. The decrease in body width may be associated with a reduced body condition previously documented in fish from sulfidic habitats [52], [87], [88].
Our morphological analyses also provided critical insights about the rediscovered P. thermalis. Most importantly, specimens collected by Heller in 1848 mostly grouped with samples from our 2012 survey, suggesting we have visited the locality described in Heller's autobiography [40]. We found P. thermalis to exhibit a typical sulfide spring body shape and to be phenotypically similar to P. sulphuraria. Qualitatively, this is also the case for color patterns (pronounced turquoise highlights on the abdomen with a relatively dark dorsal coloration, which are not found in sulfidic P. mexicana populations) and lateral lip appendages on the lower jaw (authors, personal observation), which are mentioned in the species description of P. sulphuraria [45]. Note, however, that such lip appendages are not a diagnostic trait for P. sulphuraria (s. l.), as specimens from the La Gloria population do not exhibit this morphological trait [89]. Despite the close morphological affinity of P. thermalis to P. sulphuraria (as compared to sulfidic and non-sulfidic populations of P. mexicana), our analyses indicated significant differences in body shape between the two species both in the lateral and dorsal projection.
Phylogenetic analyses and population genetics
The broad phylogenic relationships uncovered in our study match the patterns of other studies [26], [34], [65], [66]. We found P. thermalis collected in both La Esperanza springs (Ixtapangajoya drainage) to be sister to P. sulphuraria from the Baños del Azufre population (Pichucalco drainage), which together were sister to P. sulphuraria from La Gloria (also Pichucalco drainage) and formed a monophyletic group. This group (P. thermalis and P. sulphuraria together) was more closely related to the northern Mexican subspecies of P. m. limantouri than the southern P. m. mexicana populations from adjacent non-sulfidic sites, corroborating earlier investigations [34]. This suggests that colonization of sulfide springs in the Ixtapangajoya and Puyacatengo drainages by the P. m. limatouri-like ancestor shared by P. thermalis and P. sulphuraria occurred earlier than sulfide spring colonization by P. m. mexicana in the other drainages. This is reflected in higher genetic divergences between sulfide spring and adjacent non-sulfidic populations in the Ixtapangajoya and Pichucalco drainages (2.0–2.2% in mitochondrial genes) compared to the Puyacatengo and Tacotalpa drainages (0.1–0.4%; see [34] for a discussion).
Our results indicate that sulfide springs in the Pichucalco and the Ixtapangajoya drainages were not colonized independently, but rather P. sulphuraria and P. thermalis are of a single evolutionary origin despite their current distribution in independent drainages, which in the area of the sulfide springs are separated by mountainous terrain reaching more than 500 meters above the surrounding elevation. This can be explained by the dynamic nature of the courses of major river systems in southern Mexico [90], [91]. Historically, the Grijalva River was an independent deltaic system that followed the course of Ixtapangajoya river [92], presenting an opportunity for connections between currently independent tributaries, particularly during periods of heavy rain and flooding between tributaries. However, considering the reduced viability of sulfide adapted fish in non-sulfidic environments [43], [93] and consequently the absence of sulfide-adapted ecotypes even in proximate freshwater habitats [44], it remains unclear how colonization through stretches of unsuitable habitats was possible even in the presence of potential connections among drainages. Hence, the alternative hypothesis is that colonization of different springs in the two drainages could have occurred independently by a once widespread ancestor (a lineage with close affinities to extant northern Mexican P. m. limantouri) with standing genetic variation for traits adaptive to sulfidic springs. Such a scenario was recently supported in stickleback, where low rates of gene flow from freshwater to marine populations maintain freshwater alleles in the marine environments at low frequency, such that selection upon colonization of a new freshwater system can rapidly reassemble freshwater ecotypes based on allelic variants already present in the ancestral population [94], [95]. The currently available data does not allow for rigorously testing these contrasting hypotheses, and additional research including a more thorough analysis of the northern Mexican P. m. limantouri is required to elucidate historical patterns of sulfide spring colonization in the Ixtapangajoya and Puyacatengo drainages. Nonetheless, our data indicate that sulfide spring colonization may not have occurred independently in different drainages, adding an additional layer of complexity in the analysis of speciation patterns in sulfidic spring fishes.
Our population genetic analyses largely supported the phylogeny in that the uppermost level of population differentiation included two clusters distinguishing between populations in sulfidic (P. sulphuraria and P. thermalis) and non-sulfidic (P. mexicana) environments irrespective of the drainage of origin. As such, the results generally support previously uncovered patterns of strong reproductive isolation between sulfide spring residents and fish from adjacent non-sulfidic sections of the same drainage [43]. Reproductive isolation among ecotypes is at least partially mediated by natural and sexual selection against immigrants, where migrant individuals from the opposite habitat type have reduced survivability and are discriminated against during mate choice [43], [46], [94]. Our analyses also found strong support for genetic structure with K = 7 divergent clusters. At this finer scale, the two P. thermalis populations from the Ixtapangajoya drainage were clearly distinct from the two P. sulphuraria populations in the adjacent Pichucalco drainage, reflecting the absence of gene flow due to the lack of contemporary connections between the two drainages.
Taxonomic considerations and conclusions
The taxonomic history of Poecilia thermalis remained uncertain since its original description by Steindachner in 1863. Our study revisited the status of P. thermalis based on recently collected material from the type locality and museum specimens using morphological, phylogenetic, and population genetic approaches. Based on our findings, we can clearly reject the previously prevalent notion that P. thermalis Steindachner 1863 is synonymous to either P. salvatoris [38], P. sphenops [96], or P. mexicana (sensu stricto, [39]). However, in relation to the sulfide spring populations from the Pichucalco drainage, currently denominated P. sulphuraria Alvarez 1947, taxonomic change can proceed in the form of two alternatives. Sulfide spring populations from the Ixtapangajoya (P. thermalis) and the Pichucalco (P. sulphuraria) could be considered as derivatives from the same evolutionary lineage and therefore considered the same species. In this case, the older names takes precedence [97] and P. sulphuraria [45] would be designated as a junior synonym of P. thermalis. Alternatively, P. thermalis could be designated as a valid, distinct species restricted to the Ixtapangajoya drainage, which would require the name P. sulphuraria to be restricted to the type locality (Baños del Azufre) and the sulfide spring population at La Gloria (currently a population of P. sulphuraria) to be considered a distinct species awaiting formal description. This interpretation is supported by reciprocal monophyly, significant population genetic differentiation as evident from FST values and principal components analysis, and significant differences in body shape among all three groups. Examination of additional characters, especially meristic traits and the structure of the male copulatory organ (gonopodium), commonly used in poeciliid systematics will hopefully lead to the resolution of the taxonomic conundrum surrounding P. thermalis.
Regardless of the taxonomic conclusions, our study has direct implications for the conservation of the sulfide spring populations in the Ixtapangajoya and Pichucalco drainages. Currently, P. sulphuraria is listed as threatened and federally protected by the Mexican government [98]. In addition, the IUCN has listed the species as critically endangered because of a very limited distribution [97], and the species is threatened by deforestation, farming, recreational activities, and more recently by extensive palm oil plantations [89]. Despite these concerns, no conservation measures have been implemented to mitigate these effects [89]. Potential taxonomic changes will require according changes in the list of endangered species in Mexico. Whether all three populations will be designated as P. thermalis or as three distinct species in the future, they clearly represent unique evolutionary lineages with highly restricted distributions meriting separate management and a high priority for conservation [99].
Overall, this study confirms the role of hydrogen sulfide in shaping convergent, phenotypic evolution in sulfide spring fishes and causing reproductive isolation between populations residing in proximate sulfidic and non-sulfidic environments. It also illustrates how an integrative, mechanistic approach to studying phenotypic evolution and speciation can inform taxonomy.
Acknowledgments
We thank C. Passow, D. Alba, N. Franssen, and Z. Culumber for their assistance in the field, and D. Miller for his pictorial depiction of Mexican mollies. We also thank A.A. Gonzalez-Diaz for providing samples of P. butleri, as well as E. Mikschi and H. Wellendorf for providing access to the type specimens of P. thermalis. H. Lopez-Fernandez, N. Lujan, and D. Taphorn provided insightful comments about the taxonomic situation of P. thermalis. This is publication number 1459 of the Biodiversity Research and Teaching Collections and 231 of the Center for Biosystematics and Biodiversity, both at Texas A&M University.
Author Contributions
Conceived and designed the experiments: M. Palacios LAR GV MT. Performed the experiments: M. Palacios M. Plath CE HL GV MT. Analyzed the data: M. Palacios M. Plath CE HL GV MT. Contributed reagents/materials/analysis tools: LAR M. Plath AL GV MT. Wrote the paper: M. Palacios GV MT.
References
- 1. Kawecki TJ, Ebert D (2004) Conceptual issues in local adaptation. Ecol Lett 7: 1225–1241.
- 2.
Schluter D (2000) The ecology of adaptive radiation. New York, USA: Oxford University Press.
- 3. Schluter D (2001) Ecology and the origin of species. Trends Ecol Evol 16: 372–380.
- 4. Rundle HD, Nosil P (2005) Ecological speciation. Ecol Lett 8: 336–352.
- 5. Grant BR, Grant PR (1989) Natural selection in a population of Darwin's finches. Am Nat 133: 377–393.
- 6. Schluter D, McPhail JD (1992) Ecological character displacement and speciation in sticklebacks. Am Nat 140: 85–108.
- 7. Losos JB (1992) The evolution of convergent structure in Caribbean Anolis communities. Syst Biol 41: 403–420.
- 8. Lu G, Bernatchez L (1999) Correlated trophic specialization and genetic divergence in sympatric lake whitefish ecotypes (Coregonus clupeaformis): support for the ecological speciation hypothesis. Evolution 53: 1491–1505.
- 9. Vamosi SM, Schluter D (2004) Character shifts in the defensive armor of sympatric sticklebacks. Evolution 58: 376–385.
- 10. Langerhans RB, Layman CA, Shokrollahi A, DeWitt TJ (2004) Predator-driven phenotypic diversification in Gambusia affinis. Evolution 58: 2305–2318.
- 11. Magurran AE, Magurran AE (1998) Population differentiation without speciation. Philos Trans R Soc Lond B Biol Sci 353: 275–286.
- 12. Hendry AP, Nosil P, Rieseberg LH (2007) The speed of ecological speciation. Funct Ecol 21: 455–464.
- 13. Nosil P (2008) Speciation with gene flow could be common. Mol Ecol 17: 2103–2106.
- 14. Hey J, Waples RS, Arnold ML, Butlin RK, Harrison RG (2003) Understanding and confronting species uncertainty in biology and conservation. Trends Ecol Evol 18: 597–603.
- 15.
Darwin C (1859) On the origin of species by means of natural selection. London, UK: Murray.
- 16. Padial JM, Miralles A, De la Riva I, Vences M (2010) Review: The integrative future of taxonomy. Front Zool 7: 1–14.
- 17.
Dobzhansky TG, Dobzhansky T (1937) Genetics and the origin of species. New York, USA: Columbia University Press.
- 18.
Mayr E (1942) Systematics and the origin of species, from the viewpoint of a zoologist. New York, USA: Colombia University Press.
- 19. Hausdorf B (2011) Progress toward a general species concept. Evolution 65: 923–931.
- 20. Gaston KJ, Mound LA (1993) Taxonomy, hypothesis testing and the biodiversity crisis. Proc R Soc Lond B Biol Sci 251: 139–142.
- 21. Pfennig DW, Wund MA, Snell-Rood EC, Cruickshank T, Schlichting CD, et al. (2010) Phenotypic plasticity's impacts on diversification and speciation. Trends Ecol Evol 25: 459–467.
- 22.
Olvido AE, Mousseau TA (2012) Geographical Variation. eLS. DOI: https://doi.org/10.1002/9780470015902.a0001736.pub3.
- 23. Rosen DE, Bailey RM (1963) The poeciliid fishes (Cyprinodontiformes), their structure, zoogeography, and systematics. Bull Am Mus Nat Hist 126: 1–176.
- 24. Poeser FN (2011) A New Species of Poecilia from Honduras (Teleostei: Poeciliidae). Copeia 2011: 418–422.
- 25.
Reis RE, Kullander, Sven O, Ferraris, Carl J (2003) Checklist of freshwater fishes of South and Central America. Porto Alegre, Brazil: EDICPURS.
- 26. Alda F, Reina RG, Doadrio I, Bermingham E (2013) Phylogeny and biogeography of the Poecilia sphenops species complex (Actinopterygii, Poeciliidae) in Central America. Mol Phylogenet Evol 66: 1011–1026.
- 27. Poeser FN (2003) Geographic variation in Poecilia Bloch and Schneider, 1801 (Teleostei: Poeciliidae), with descriptions of three new species and designation of lectotypes for P. dovii Günther, 1866 and for P. vandepolli van Lidth de Jeude, 1887. Proc Biol Soc Wash 116: 356–379.
- 28. Hubbs CL (1926) Studies of the fishes of the order cyprinodontes VI. University of Michigan, Museum of Zoology, Miscellaneous Publications 16: 1–86.
- 29. Regan CT (1913) A revision of the cyprinodont fishes of the subfamily Poeciliinae. Proc Zool Soc Lond 2: 977–1018.
- 30.
Seda JB (2010) Intra- and interspecific variation in morphology and behaviors in the sailfin molly, Poecilia latipinna (Poeciliidae: Poecilia: Mollienesia). Master's thesis. Clemson University, Clemson, South Carolina.
- 31. Zúñiga-Vega JJ, Suárez-Rodríguez M, Espinosa-Pérez H, Johnson JB (2011) Morphological and reproductive variation among populations of the Pacific molly Poecilia butleri. J Fish Biol 79: 1029–1046.
- 32. Parzefall J (1989) Sexual and Aggressive Behaviour in Species Hybrids of Poecilia mexicana and Poecilia velifera (Pisces, Poeciliidae). Ethology 82: 101–115.
- 33. Kittell MM, Harvey MN, Contreras-Balderas S, Ptacek MB (2005) Wild-caught hybrids between sailfin and shortfin mollies (Poeciliidae, Poecilia): morphological and molecular verification. Hidrobiologica 15: 131–137.
- 34. Tobler M, Palacios M, Chapman LJ, Mitrofanov I, Bierbach D, et al. (2011) Evolution in extreme environments: replicated phenotypic differentiation in livebearing fish inhabiting sulfidic springs. Evolution 65: 2213–2228.
- 35. Steindachner F (1863) Beiträge zur Kenntniss der Sciaenoiden Brasiliens und der Cyprinodonten Mejicos. Sitzungsberichte der Mathematisch-Naturwissenschaftlichen Classe der Kaiserlichen Akademie der Wissenschaften 48: 162–185.
- 36.
Günther A (1866) Catalogue of the Physostomi, containing the families Salmonidae, Percopsidae, Galaxidae, Mormyridae, Gymnarchidae, Esocidae, Umbridae, Scombresocidae, Cyprinodontidae, in the collection of the British Museum. Catalogue of fishes in the British Museum. London, UK: Taylor and Francis.
- 37.
Poeser F (2003) From the Amazon river to the Amazon molly and back again. PhD dissertation. Amsterdam University, Amsterdam, NL: IBED.
- 38. Regan CT (1907) Diagnoses of new Central-American freshwater fishes of the families Cyprinodontidae and Mugilidae. Ann mag nat hist 19: 64–66.
- 39. Miller RR (1983) Checklist and Key to the Mollies of Mexico (Pisces: Poeciliidae: Poecilia, Subgenus Mollienesia). Copeia 1983: 817–822.
- 40.
Heller KB (1853) Reisen in Mexico in den jahren 1845-1848. Leipzig, Germany: Verlag Von Wilhelm Engelmann.
- 41. Plath M, Tobler M, Riesch R, García de León FJ, Giere O, et al. (2007) Survival in an extreme habitat: the roles of behaviour and energy limitation. Naturwissenschaften 94: 991–996.
- 42.
Pfenninger M, Lerp H, Tobler M, Passow C, Kelley JL, et al.. (2009) Parallel evolution of cox-genes in H2S-tolerant fish as key adaptation to a toxic environment. Submitted.
- 43.
Plath M, Pfenninger M, Lerp H, Riesch R, Eschenbrenner C, et al.. (2013) Genetic differentiation and selection against migrants in evolutionarily replicated extreme environments. Evolution in press.
- 44. Plath M, Hermann B, Schröder C, Riesch R, Tobler M, et al. (2010) Locally adapted fish populations maintain small-scale genetic differentiation despite perturbation by a catastrophic flood event. BMC Evol Biol 10: 256.
- 45. Alvarez del Villar J (1947) Descripción de una nueva especie de Mollienisia capturada en Baños del Azufre, Tabasco (Pisces, Poeciliidae). An Esc Nac Cienc Biol Mex 5: 275–281.
- 46. Plath M, Hauswaldt JS, Moll K, Tobler M, García de León FJ, et al. (2006) Local adaptation and pronounced genetic differentiation in an extremophile fish, Poecilia mexicana, inhabiting a Mexican cave with toxic hydrogen sulphide. Mol Ecol 16: 967–976.
- 47. Tobler M, Riesch R, Tobler CM, Schulz-Mirbach T, Plath M (2009) Natural and sexual selection against immigrants maintains differentiation among micro-allopatric populations. J Evol Biol 22: 2298–2304.
- 48. Tobler M, DeWitt TJ, Schlupp I, García de León FJ, Herrmann R, et al. (2008) Toxic hydrogen sulfide and dark caves: phenotypic and genetic divergence across two abiotic environmental gradients in Poecilia mexicana. Evolution 62: 2643–2659.
- 49. Armienta MA, De la Cruz-Reyna S (1995) Some hydro-geochemical fluctuations observed in Mexico related to volcanic activity. Appl Geochem 10: 215–227.
- 50. Rosales Lagarde L, Boston PJ, Campbell A, Stafford KW (2006) Possible structural connection between Chichón Volcano and the sulfur-rich springs of Villa Luz Cave (aka Cueva de las Sardinas), Southern Mexico. Association for Mexican Cave Studies Bulletin 19: 177–184.
- 51.
Rosales Lagarde L (2012) Investigations of karst brackish-sulfidic springs and their role in sulfidic springs and their hydrogeology, subsurface water-rock interactions, and speleogenesis at the northern Sierra de Chiapas, Mexico. PhD dissertation. New Mexico Institute of Mining and Technology, Socorro, New Mexico.
- 52. Tobler M, Schlupp I, Heubel KU, Riesch R, de León FJG, et al. (2006) Life on the edge: hydrogen sulfide and the fish communities of a Mexican cave and surrounding waters. Extremophiles 10: 577–585.
- 53. Roach KA, Tobler M, Winemiller KO (2011) Hydrogen sulfide, bacteria, and fish: a unique, subterranean food chain. Ecology 92: 2056–2062.
- 54. Tobler M, Hastings L (2011) Convergent patterns of body shape differentiation in four different clades of poeciliid fishes inhabiting sulfide springs. Evol Biol 38: 412–421.
- 55.
Rohlf FJ (2010) tpsDig2 version 2.16. Computer program and documentation distributed through http://life.bio.sunysb.edu/morph/. Accessed 2012 July 8.
- 56.
Zelditch ML, Swiderski DL, Sheets HD (2012) Geometric morphometrics for biologists: a primer. New York, USA: Academic Press.
- 57.
Rohlf FJ (2007) tpsRelw version 1.45. Computer program and documentation distributed through http://life.bio.sunysb.edu/morph/. Accessed 2012 July 8.
- 58. Langerhans RB, DeWitt TJ (2004) Shared and unique features of evolutionary diversification. Am Nat 164: 335–349.
- 59. Langerhans RB (2009) Trade-off between steady and unsteady swimming underlies predator-driven divergence in Gambusia affinis. J Evol Biol 22: 1057–1075.
- 60.
Rohlf FJ (2005) tpsRegr version 1.38. Computer program and documentation distributed through http://life.bio.sunysb.edu/morph/. Accessed 2012 July 8.
- 61. Saitou N, Nei M (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol 4: 406–425.
- 62.
Schmidt TR, Bielawski JP, Gold JR (1998) Molecular phylogenetics and evolution of the cytochrome b gene in the cyprinid genus Lythrurus (Actinopterygii: Cypriniformes). Copeia: 14–22.
- 63.
Broughton RE, Gold JR (2000) Phylogenetic relationships in the North American cyprinid genus Cyprinella (Actinopterygii: Cyprinidae) based on sequences of the mitochondrial ND2 and ND4L genes. Copeia: 1–10.
- 64. Kocher TD, Conroy JA, McKaye KR, Stauffer JR, Lockwood SF (1995) Evolution of NADH dehydrogenase subunit 2 in East African cichlid fish. Mol Phylogenet Evol 4: 420–432.
- 65. Meredith RW, Pires MN, Reznick DN, Springer M (2010) Molecular phylogenetic relationships and the evolution of the placenta in Poecilia (Micropoecilia) (Poeciliidae: Cyprinodontiformes). Mol Phylogenet Evol 55: 631–639.
- 66. Meredith RW, Pires MN, Reznick DN, Springer MS (2011) Molecular phylogenetic relationships and the coevolution of placentotrophy and superfetation in Poecilia (Poeciliidae: Cyprinodontiformes). Mol Phylogenet Evol 59: 148–157.
- 67. Katoh K, Toh H (2008) Recent developments in the MAFFT multiple sequence alignment program. Brief. Bioinform. 9: 286–298.
- 68.
Nylander JAA (2004) MrModeltest v2. Program distributed by the author. Evolutionary Biology Centre, Uppsala University. Available from http://www.abc.se/~nylander/mrmodeltest2/mrmodeltest2.html. Accessed 2013 January 8.
- 69. Stamatakis A (2006) RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22: 2688–2690.
- 70. Stamatakis A, Hoover P, Rougemont J (2008) A rapid bootstrap algorithm for the RAxML web servers. Syst Biol 57: 758–771.
- 71.
Zwickl DJ (2006) Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasers under the maximum likelihood criterion. PhD dissertation. The University of Texas at Austin, Austin, Texas.
- 72. Sukumaran J, Holder MT (2010) DendroPy: a Python library for phylogenetic computing. Bioinformatics 26: 1569–1571.
- 73. Ronquist F, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.
- 74. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, et al. (2012) MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol 61: 539–542.
- 75.
Rambaut A, Drummond AJ (2007) Tracer version 1.4. Computer program and documentation distributed by the author. Available from http://beast.bio.ed.ac.uk/Tracer. Accessed 2013 January 8.
- 76. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, et al. (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28: 2731–2739.
- 77. Tiedemann R, Moll K, Paulus KB, Schlupp I (2005) New microsatellite loci confirm hybrid origin, parthenogenetic inheritance, and mitotic gene conversion in the gynogenetic Amazon molly (Poecilia formosa). Mol Ecol Notes 5: 586–589.
- 78.
Slattery P, Eschenbrenner C, Arias-Rodriguez L, Streit B, Bierbach D, et al.. (2012) Twelve new microsatellite loci for the sulphur molly (Poecilia sulphuraria) and the related Atlantic molly (P. mexicana). Conserv Genet Resour: 1–3.
- 79. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959.
- 80. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14: 2611–2620.
- 81.
Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour: 1–3.
- 82. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes 6: 288–295.
- 83. Tobler M, Riesch R, Tobler CM, Plath M (2009) Compensatory behavior in response to sulfide-induced hypoxia affects time budgets, feeding efficiency, and predation risk. Evol Ecol Res 11: 935–948.
- 84. Winemiller KO (1989) Development of dermal lip protuberances for aquatic surface respiration in South American characid fishes. Copeia 1989: 382–390.
- 85. Scarabotti PA, Julieta Parma M, López JA, Ghirardi R (2009) Dermal lip protuberances associated with aquatic surface respiration in juveniles of the piscivorous characid Salminus brasiliensis (Actinopterygii: Characidae). Neotrop Ichthyol 7: 459–464.
- 86. Scarabotti PA, López JA, Ghirardi R, Parma MJ (2011) Morphological plasticity associated with environmental hypoxia in characiform fishes from neotropical floodplain lakes. Environ Biol Fishes 92: 391–402.
- 87. Plath M, Heubel KU, García de León FJ, Schlupp I (2005) Cave molly females (Poecilia mexicana, Poeciliidae, Teleostei) like well-fed males. Behav Ecol Sociobiol 58: 144–151.
- 88. Tobler M (2008) Divergence in trophic ecology characterizes colonization of extreme habitats. Biol J Linn Soc 95: 517–528.
- 89. Tobler M, Plath M (2009) Threatened fishes of the world: Poecilia sulphuraria (Alvarez, 1948) (Poeciliidae). Environ Biol Fishes 85: 333–334.
- 90. Psuty NP (1965) Beach-ridge development in Tabasco, Mexico. Ann Assoc Am Geog 55: 112–124.
- 91.
West RC, Psuty NP, Thom BG (1969) The Tabasco lowlands of southeastern Mexico. Louisiana, USA: Louisiana State University Press.
- 92.
Böse E (1905) Reseña acerca de la geología de Chiapas y Tabasco (Vol. 20). Mexico: Imprenta y fototipía de la Secretaría de Fomento.
- 93. Plath M, Riesch R, Oranth A, Dzienko J, Karau N, et al. (2010) Complementary effect of natural and sexual selection against immigrants maintains differentiation between locally adapted fish. Naturwissenschaften 97: 769–774.
- 94. Schluter D, Conte GL (2009) Genetics and ecological speciation. Proc Natl Acad Sci 106: 9955–9962.
- 95. Jones FC, Grabherr MG, Chan YF, Russell P, Mauceli E, et al. (2012) The genomic basis of adaptive evolution in threespine sticklebacks. Nature 484: 55–61.
- 96. Gordon MS, Rosen DE (1962) A cavernicolous form of the poeciliid fish Poecilia sphenops from Tabasco, Mexico. Copeia 1962: 360–368.
- 97.
IUCN (2012) IUCN Red List of Threatened Species. Version 2012.2. www.iucnredlist.org. Accesed 2013 Jan 15.
- 98.
Sedesol (2010) NORMA Oficial Mexicana NOM-059-SEMARNAT-2010, Protección ambiental-Especies nativas de México de flora y fauna silvestres-Categorías de riesgo y especificaciones para su inclusión, exclusión o cambio-Lista de especies en riesgo. Diario Oficial de la Federacion 2: 1–78. http://dof.gob.mx/nota_detalle.php?codigo=5173091&fecha=30/12/2010 Accesed 2013 Jan 15.
- 99. Ryder OA (1986) Species conservation and systematics: the dilemma of subspecies. Trends Ecol Evol 1: 9–10.