Population level genetic divergence and phylogenetic placement of Mexican shortfin mollies (Mollienesia: Poecilia: Poeciliidae)

Abstract Mexico is a megadiverse region with a complex geological history, but it remains unclear to what extent the distribution of freshwater fish has been influenced by geographic barriers. This study examines the population level genetic divergence and phylogenetic relationships of species in the shortfin group of the subgenus Mollienesia (genus Poecilia), a group of live-bearing fishes that are widely distributed across Mexico, with sampling at a small geographic scale. Samples from over 50 locations were analyzed for six species by using phylogenetic and haplotype network approaches to assess genetic diversity across geographic ranges and to refine the distributions of species in this group. The results indicate that Mexican species have diversified following multiple, independent invasions from Middle America. Two species found north of the Trans-Mexican Volcanic Belt (TMVB) and one transversal species exhibited weak phylogenetic structure, likely due to the lack of physiographic barriers, recent colonization, and high dispersal rates among regions. In contrast, three species found south of the TMVB exhibited strong phylogenetic structure, reflecting a longer presence in the area and multiple physiographic barriers that isolated populations. This study identified mechanisms driving divergence and speciation, expanded the known range of several species, and resolved taxonomic uncertainties of populations.


INTRODUCTION
Mexico's complex geological history has led to diverse ecosystems, making the region rich in flora and fauna.The historical geological activity is reflected in a landscape dominated by mountain chains that give rise to isolated hydrological units and regions with starkly different ecological conditions (Domínguez-Domínguez et al., 2006).In the north of the country, the mountain ranges of the Sierra Madre Oriental and the Sierra Madre Occidental extend along the eastern and western coast, respectively, and in the south, they join the Trans-Mexican Volcanic Belt (TMVB; Fig. 1), which transects the country at a latitude of about 20 ºN and marks the southern edge of the North American plate (Ferrari et al., 1999).In the southern region of Mexico, the Sierra Madre del Sur and the Sierra Madre de Chiapas are separated by the Isthmus of Tehuantepec (Fig. 1).These mountain ranges produce a wide range of geographical landforms, including high mountains with permanent snow cover, highland plateaus, deep valleys, and coastal plains (Contreras-Balderas et al., 2008).Mexico is also unique in encompassing organisms from both the Neartic (North America south to the Mexican plateau) and Neotropical (Mexican central coasts to South America) biogeographic zones (Morrone, Márquez, 2001).The combination of the geological diversity and the interface of different biogeographic zones make Mexico a megadiverse region in terms of ecological habitats, species diversity, and levels of endemism (Conabio, 1998;Metcalfe et al., 2000;Hufnagel, Mics, 2021).Mexican molly phylogeography the subgenus Mollienesia (P.butleri, P. limantouri Jordan & Snyder, 1899, P. mexicana, P. nelsoni, P. sphenops, P. sulphuraria, and P. thermalis) from 66 locations across Mexico to (1) determine the number of colonization events of Mexico from nuclear Central America, (2) assess patterns of phylogeography and haplotype diversity in different species, and (3) contribute to a better understanding of the distribution of species in this group.

MATERIAL AND METHODS
Specimen acquisition.We sampled specimens of the subgenus Mollienesia throughout the subgenus's distribution in Mexico (Fig. 1; see Tab.S1 for locality data).Fish were captured using electrofishing, seines, cast and dip nets.Immediately after capture, some sampled fish were euthanized with buffered MS222; others were sampled and released (per permit restrictions).Regardless of sampling method, fin clips (right pectoral fin) were preserved in 95% ethanol for molecular analyses, and specimens were fixed in a 10% formaldehyde solution following protocols approved by the Texas A&M University and Oklahoma State University Committees on Use and Care of Animals (IACUC 2011-118 & ACUP AS10-15).Ethanol preserved tissues, DNA extractions, and formalin fixed specimens are housed in the Biodiversity, Research and Teaching Collections of the Department of Wildlife & Fisheries Sciences at Texas A&M University, College Station, Texas, USA, and the El Colegio de la Frontera Sur, San Cristóbal, Chiapas, Mexico.Voucher numbers for new specimens are provided in Tab.S1, where possible; many of our samples were fin clips (per permit restrictions) and some were used for life history dissections, and as such were not suitable as voucher specimens.Samples collected for this study were supplemented with sequence data obtained from GenBank to examine the relationship of Mollienesia species in Mexico to other members of the genus from North and Central America.We included sequence data for multiple genes (see below) for 18 out of 26 recognized species in the subgenus Mollienesia (Mateos, 2005;Alda et al., 2013;Ho et al., 2016;Palacios et al., 2016): P. kykesis, P. latipinna, P. latipunctata Meek, 1904, and P. velifera (Regan, 1914) (members of the sailfin molly clade); P. catemaconis, P. chica, P. marcellinoi, and P. sphenops (members of the P. sphenops species complex of the shortfin molly clade), as well as P. butleri, P. gillii (Kner, 1863), P. hondurensis Poeser, 2011, P. mexicana, P. nelsoni, P. orri Fowler, 1943, P. petenensis Günther, 1866 (sometimes referred to as P. gracilis), P. salvatoris Regan, 1907, P. sulphuraria, and P. thermalis (members of the P. mexicana species complex of the shortfin molly clade).In addition, we supplemented the cytochrome b dataset generated for our focal species with additional sequences.This expanded the distributional coverage for P. butleri and P. nelsoni (Zúñiga-Vega et al., 2014), spanning their entire ranges along the Pacific coast of Mexico.The distribution coverage was also extended north of our sampling for P. limantouri (Stöck et al., 2010) and south of our sampling in P. mexicana (Tobler et al., 2011).
Alignment and phylogenetic analyses.The data were partitioned into three datasets: (1) a concatenated dataset (cyt b, ND2, and RAG1, 3748 base pairs), (2) the nuclear gene alone (RAG1 1561 bp), and (3) cyt b alone (1140 bp).A strictly mitochondrial dataset (cyt b and ND2) was not analyzed because samples missing sequences for ND2 reduced the size of the dataset.All datasets were analyzed separately using FaBox online (Villesen, 2007) to determine the total number of haplotypes.The cyt b dataset was reduced from 953 sequences to 475 (139 haplotypes for the 6 species of interest).We used Partition Finder (Lanfear et al., 2012) to determine the best partition scheme and the most likely model of DNA substitution among 24 candidate models on a fixed BioNJ-JC tree based on the Bayesian information criterion (BIC; Tab.S2), separately for each gene.Phylogenetic analyses for the concatenated and the cyt b datasets were conducted with a maximum likelihood (ML) approach as implemented in RAxML GUI v. 1.0 (Stamatakis, 2006, 2008), with 500 Rapid Bootstrap searches followed by an ML search.The complex general time reversible (GTR) + Γ (gamma distribution for rate variation among sites) model was chosen for each partition.The bootstrap trees were summarized with a Sumtrees script that used a 50% percent majority rule consensus parameter in DendroPy v. 3.10.1 (Sukumaran, Holder, 2010), and the final tree was rooted and visualized in FigTree v. 1.4.2(http://tree.bio.ed.ac.uk/software/figtree/).Bayesian analyses were conducted in MrBayes v. 3.2.1 (Ronquist, Huelsenbeck, 2003;Ronquist et al., 2012), implementing two runs with four chains under default parameters for 50 million generations, sampling every 10 generations.A 25% burn-in was applied and stable posterior probability values examined in Tracer v. 1.5 (Rambaut et al., 2018).Pairwise genetic distances were calculated under the Kimura-2 parameter model in MEGA v. 7 (Kumar et al., 2016) with pairwise deletion for missing data.
A haplotype network was generated with a statistical parsimony analysis of the cyt b gene dataset for the six focal species in TCS v. 1.2.1 (Clement et al., 2000), as implemented in the Popart program (Leigh, Bryant, 2015).This approach calculated the number of significant substitutions connecting haplotypes in a network by applying the algorithm developed by Templeton et al. (1992).A connection limit was set to the default of 95% to generate the haplotype networks between closely related sequences.

Phylogenetic relationships of Mexican mollies to other species in the subgenus
Mollienesia.To determine the number of colonization events from nuclear Central America into Mexico, we generated a phylogenetic tree with 18 of the 26 recognized species in the subgenus Mollienesia, and overall, there was high concordance between maximum-likelihood and Bayesian analyses.The phylogeny of the concatenated dataset showed moderate to high support for the two clades representing the P. sphenops and P. mexicana species complexes (Fig. 2).Most of the phylogenetic signal came from the mitochondrial genes.Analysis of RAG1 alone resulted in poorly supported phylogenetic relationships due to the low levels of sequence variation in our samples (1-2 base pairs across the entire gene; results not shown).The genetic distances across the two complexes ranged from 6.9-10.5%,while they ranged from 1.0-8.2%within the P. sphenops species complex and 0.9-6.0%within the P. mexicana complex (Tab.S3).
The P. sphenops species complex was strongly (Bayesian, > 95%) or moderately (maximum-likelihood, 74%) supported and included three species from Mexico (P.sphenops and the two endemic species P. chica and P. catemaconis), which were interspersed among three species from Middle America (P.marcellinoi and two putatively undescribed species, P. "sp.Tipitapa" and P. "sphenops sp1"; Bagley et al., 2015) (Fig. 2).Within P. sphenops, our phylogeny did not provide support for the relationships among the sampled locations that cover both the Atlantic and Pacific versants of Mexico.
The P. mexicana species complex was strongly supported (>95% with both methods) and included five Mexican and five Middle American species (Fig. 2).Poecilia hondurensis and P. orri from Middle America represented the most basal divergences in the complex (Fig. 2).The next divergence in the complex was strongly supported and included two species from Mexico's Pacific slope, P. butleri and P. nelsoni.Next, our results recovered a strongly supported P. limantouri clade on the Atlantic slope of Mexico (which lacks resolution among different sampling sites) as sister to the sulfide spring endemics P. sulphuraria and P. thermalis (Fig. 2).The P. limantouri/P.sulphuraria-P.thermalis clade is sister to a clade comprised of two putatively undescribed lineages from Middle America (P."sp 1" and P. "sp 2"; Palacios et al., 2016), P. salvatoris (also from Middle America), and the strongly/ moderately supported P. mexicana clade.Poecilia mexicana is split into three clades.Two of these clades are of samples from Mexico, with one strongly supported clade comprised of the most northern sampling locations, from the Trans-Mexican Volcanic Belt north to the state of Nuevo León.The other clade includes the remaining samples, which come from south of the Trans-Mexican Volcanic Belt and span both sides of the Isthmus of Tehuantepec (Fig. 2).The two Mexican clades are not sisters, as a lineage consisting of samples from Honduras is more closely related to the southern clade.

Phylogeography of Poecilia sphenops.
To assess patterns of phylogeography and haplotype diversity within different species, we generated haplotype networks and observed potential patterns of differentiation across known geographic barriers.The widespread species P. sphenops is found across the Atlantic and Pacific versants in Mexico, and genetic distances within this species ranged from 0 to 1.0% (Tab.S4).Our analyses revealed three divergent groups.The first of these (Texistepec) is represented by a single haplotype from the Rio Coatzacoalcos of the southern Atlantic coast (state  Maura Palacios, Alfonso A. González-Díaz, Lenin A. Rodriguez, Mariana Mateos, Rocío Rodiles-Hernández, Michael Tobler, Gary Voelker of Veracruz), just west of the Isthmus of Tehuantepec (Fig. 3).The second group (clade A1) was composed of one Pacific slope haplotype from the Guerrero coast west of the Isthmus of Tehuantepec, and two haplotypes from the Atlantic versant in Chiapas west of the Isthmus of Tehuantepec (Fig. 3).The main P. sphenops clade (clade A2 in Fig. 3) included one widely distributed haplotype, found along the Atlantic coast north of the Trans-Mexican Volcanic Belt, in the Balsas depression, and along the Pacific coast east and west of the Isthmus of Tehuantepec.Clade A2 also exhibited three strongly supported sub-clades, the first consisting of one haplotype from the Balsas depression and one from a Pacific location at the Isthmus of Tehuantepec (Rio Ostuta on the Oaxacan coast).The second sub-clade consisted of four Pacific haplotypes from the Chiapas coast east of the Isthmus of Tehuantepec, and the third sub-clade was comprised of three haplotypes from the Atlantic versant east of the Isthmus of Tehuantepec in Chiapas (Fig. 3).Additionally, we found five other P. sphenops haplotypes with unresolved relationships; these samples were geographically dispersed across the Balsas depression and the Pacific versant on both sides of the Isthmus of Tehuantepec (Fig. 3).
The haplotype network of P. sphenops indicated that haplotypes were separated between 1-10 mutational steps.The most common haplotype included samples from the Pacific coast, including both sides of the Isthmus of Tehuantepec and the Balsas depression, and one sample from Atlantic versant north of the Trans-Mexican Volcanic Belt.The most divergent haplotypes in P. sphenops were found along the Pacific coast west of the Isthmus of Tehuantepec and in rivers of the Atlantic versant, both west (Rio Coatzocoalcos) and east (Chiapas highlands) of the Isthmus (Fig. 3).The remaining haplotypes, with two exceptions, were all found on the Pacific versant, including the Balsas depression.Overall, the lack of clear phylogeographic patterns and haplotype sharing across major biogeographic boundaries indicate low levels of differentiation within P. sphenops.
Phylogeography of Poecilia butleri and P. nelsoni.The sister taxa P. butleri and P. nelsoni are found on either side of the Trans-Mexican Volcanic Belt.The genetic distances within both P. butleri and P. nelsoni ranged from 0-2.0%, and the two species differed by 4.6% (Tab.S4).
Poecilia butleri was restricted to the Pacific versant north of the Trans-Mexican Volcanic Belt.The phylogenetic structure for this species showed an initial divergence represented by a single haplotype from our most northern sample in Sinaloa (Rio Choix).The next divergence (clade B2 in Fig. 4) included a clade composed of two haplotypes from the southernmost portion of the P. butleri range (San Blas and Rio El Palillo in southern Nayarit).Lastly, the strongly supported clade B1 contained the remaining haplotypes, which formed a five lineage polytomy with no obvious geographic structuring (Fig. 4).Within clade B1, two lineages represented single haplotypes, one from the northern (Sinaloa) and one from the southern portion (Nayarit) of P. butleri's range.Two other lineages represented moderately to strongly supported sub-clades comprising six and three haplotypes from Sinaloa.The fifth lineage was a large weakly supported clade of individuals collected from both the northern and southern portions of the range.The haplotype network of P. butleri had a maximum of 17 mutational steps between haplotypes, with most haplotypes differing between 1-6 steps (Fig. 4).The most divergent haplotypes were from Rio Choix (192), the most northern sampling locality in Sinaloa, and from El Palillo (456), San Blas (399), and Acaponeta (408), all southern localities in Nayarit.Poecilia nelsoni is distributed on the Pacific versant south of the Trans-Mexican Volcanic belt, spanning both sides of the Isthmus of Tehuantepec.Although all P. nelsoni samples were recovered as monophyletic, the basal node consisted of a polytomy of six lineages (Fig. 4).Two lineages consisting of a single haplotype, a third lineage representing two haplotypes (clade C2), and a fourth lineage representing many haplotypes (clade C3) were mostly collected west of the Isthmus of Tehuantepec in Guerrero, with one individual from the Balsas depression (Fig. 4).Of the two remaining P. nelsoni lineages, one (clade C1) was composed almost entirely of individuals collected just south of the Trans-Mexican Volcanic Belt, including sites in Jalisco, Colima, and Michoacan, and one individual from the Balsas depression.Finally, clade C4 was composed of individuals collected west of the Isthmus of Tehuantepec in Oaxaca and individuals collected south of the Isthmus of Tehuantepec in Chiapas.The haplotype network for P. nelsoni underscored the phylogenetic results, in that multiple divergent haplotype groups were found at the northern extent of P. nelsoni's range, and along the coast of Guerrero and Oaxaca.In addition, haplotypes were shared across the Isthmus of Tehuantepec, and the Balsas River samples were more closely associated with samples from other regions than with each other.Phylogeography of Poecilia mexicana sensu lato.The closely related species P. limantouri, P. thermalis and P. sulphuraria, and P. mexicana formed, along with several Middle American species, the last major clade in the phylogeny (Figs. 2, 5).Among these species, the hydrogen-sulfide spring endemics P. thermalis and P. sulphuraria were most genetically similar (0-0.5% sequence divergence within each species), whereas P. mexicana (within species range: 0-2.0%) and P. limantouri (within species range: 0-1.7%) were more genetically diverse (Tab.S4).Poecilia limantouri was predominantly found north of the TMVB, but with a small range (represented by two sampling points) just south of this geographic barrier.Most haplotypes of P. limantouri across the distributional range did not exhibit any phylogenetic clustering or lacked support from the maximum likelihood analyses.The single sub-clade with strong support (D1) included two haplotypes: one from the state of Hidalgo and the other from the southern sampling points in the state of Veracruz (Fig. 5).The haplotype network exhibited between 1-6 mutational steps between haplotypes.The species P. sulphuraria and P. thermalis from the states of Tabasco and Chiapas were recovered in two sub-clades with moderate support.The first (E1) included individuals of P. sulphuraria from the Baños del Azufre site and of P. thermalis from the La Esperanza spring complex.The second sub-clade (E2) was formed exclusively of haplotypes of P. sulphuraria from the La Gloria spring complex.The haplotype network showed almost no variation in these taxa, and interestingly the Baños del Azufre population of P. sulphuraria clustered more closely with P. thermalis (clade E1) than the La Gloria population of P. sulphuraria (clade E2; Fig. 5).
Finally, phylogenetic analyses recovered a monophyletic P. mexicana, with further phylogenetic structuring of varying support (Fig. 5).Maura Palacios, Alfonso A. González-Díaz, Lenin A. Rodriguez, Mariana Mateos, Rocío Rodiles-Hernández, Michael Tobler, Gary Voelker individuals collected just south of TMVB, and two sites just west of the Isthmus of Tehuantepec in the Los Tuxlas region (clade F1 in Fig. 5).Another sub-clade represented individuals collected from the Atlantic versant of Honduras, and this sub-clade was recovered as sister to the last sub-clade (Clade F2), comprising individuals collected from Tabasco and Chiapas, except for one individual collected in Veracruz (Fig. 5).The haplotype network showed a range of variation with 1-15 mutational steps between haplotypes.
Our analyses indicated extensions of the documented ranges for P. mexicana and P. limantouri, which were previously believed to either occur only south and only north of the TMVB, respectively (Miller et al., 2005).Both species' distributions appear to span the TMVB (Fig. 5), and the northern haplotypes in P. mexicana represent a highly divergent lineage that may represent a yet unrecognized species that also differs in a series of morphological traits from P. mexicana south of the TMVB (M.Tobler, unpublished data).We also found the southern range of Poecilia nelsoni in Mexico to extend into southern Oaxaca and Chiapas.

DISCUSSION
The aim of this study was to provide a better understanding of the patterns of divergence and speciation across the geologically complex Mexican landscape in freshwater fishes of the subgenus Mollienesia.The overall phylogenetic patterns observed among species (Fig. 2) revealed multiple colonization events from lower Central America into Mexico and strong geographic structuring by physiographic barriers and by river basins.The exception was the widely distributed P. sphenops, which was distributed across all geographic barriers discussed here (Figs. 1, 3), thereby overlapping with the distribution of other species in the subgenus Mollienesia.In contrast, P. nelsoni and P. mexicana showed signals of population structure but not in relation to the Isthmus of Tehuantepec (Figs. 4,5).This phylogeographic structure over small spatial scales, albeit weak at present, is likely the result of limited genetic exchange across smaller geologic barriers and disconnected river basins.Our results indicate that different species in the subgenus Mollienesia-despite their phenotypic similarities-vary in their vagility, and mechanistic studies are needed to understand how morphological, physiological, and behavioral differences among species may contribute to variation in dispersal.

Phylogenetic patterns of Mexican species. The phylogeny of the subgenus
Mollienesia shows an interesting pattern with respect to the Mexican species in the P. sphenops and P. mexicana species complexes.Most Mexican species in the subgenus were directly derived from Central American lineages, demonstrating the complexity of the evolution of species in this group.In the P. sphenops species complex, there were three potential invasions that led to the speciation of Mexican species (P.catemaconis, P. chica, and P. sphenops).The differentiation of P. catemaconis on the Atlantic slope of Mexico from the ancestral population was a result of the geographic isolation as Lake Catemaco became geologically separated from other hydrological units by the Sierra de los Tuxtlas (Palacios et al., 2016).The Pacific P. chica possibly also evolved due to volcanic activity isolating the Jalisco block between the Pliocene and Quaternary (Rosas-Elguera et al., 1996).

Mexican molly phylogeography
In the P. mexicana species complex, Mexican species also showed multiple independent invasions from Central America.The Pacific slope species, P. butleri and P. nelsoni, evolved in Mexico due to the uplift of the Trans Mexican Volcanic Belt, isolating populations north and south of this barrier (Mateos, 2005;Zúñiga-Vega et al., 2014).On the Atlantic versant, the evolution of species was a combination of vicariance (P.limantouri from P. petenensis in Central America; Palacios et al., 2016) and adaptation to local environmental conditions in toxic, hydrogen-sulfide-rich springs (P.sulphuraria and P. thermalis; Palacios et al., 2013).Lastly, the most recent invasion of P. mexicana may have included two independent invasion events from Central America; one represented by lineages occurring in the northern distribution and the other into the southern portion in Mexico.
Within species, strong phylogeographic structuring was observed in two species (P.nelsoni and P. mexicana) over small spatial scales.This is likely due to the fragmentation and disruption of connections between drainages, potentially by the presence of physiographic barriers such as the Balsas Depression, thereby limiting dispersal and isolating populations.This pattern is common in other Neotropical fish species distributed in Mexico, such as the genera Astyanax Baird & Girard, 1854 (Coghill et al., 2014), Profundulus Hubbs, 1924(Morcillo et al., 2016), and Poeciliopsis Regan, 1913(Mateos et al., 2002), although other species, like P. sphenops, appear to be able to cross the very same barriers.

Geological barriers influencing diversification in Mexican mollies. The Trans
Mexican Volcanic Belt (TMVB) is a physiographic feature that formed west to east between 25-2.5 million years ago (Miller et al., 2005).The TMVB is composed of ridges and volcanoes that decrease eastwardly in altitude toward Veracruz (Miller et al., 2005).This barrier is a phylogeographic break for both terrestrial and aquatic species (Huidobro et al., 2006), and our results suggest that the TMVB was associated with a vicariant speciation event on the Pacific versant, resulting in the evolution of sister species north (P.butleri) and south (P.nelsoni) of this barrier.This pattern has been previously observed in these sister taxa (Mateos, 2005;Zúñiga-Vega et al., 2014) and in other vertebrate species (Devitt, 2006;Blair, Sánchez-Ramírez, 2016;Light et al., 2016).However, the most common recent ancestors of several fish groups distributed across the TMVB originated north of the TMVB (Pérez-Rodríguez et al., 2015), whereas our study suggests an origin south of the TMVB because of the higher genetic variation in the populations of P. nelsoni (observed in most southern populations) and due to the phylogenetic structure pointing toward repeated colonization from Central America.A similar south to north crossing of the TMVB has recently been documented for Poeciliopsis scarli Meyer, Riehl, Dawes & Dibble, 1985(Conway et al., 2019).Invasion by or isolation of populations of Neotropical freshwater fish to the north of the TMVB on the Pacific Coast has also been documented in fish of the genera Astyanax (Strecker et al., 2012) andPoeciliopsis (Mateos et al., 2002).
On the Atlantic slope, the eastern most point of TMVB is called Punta del Morro (PDM) and serves as a phylogeographic break (Contreras-Balderas et al., 1996).The uplift of this barrier has caused vicariant events evident in species pairs of various families of Neartic and Neotropical freshwater fish (e.g., Lepisosteidae, Clupeidea, Cichlidae, Characidae, and Poeciliidae;Contreras-Balderas et al., 1996;Hulsey et al., 15/21 ni.bio.br| scielo.br/niMaura Palacios, Alfonso A. González-Díaz, Lenin A. Rodriguez, Mariana Mateos, Rocío Rodiles-Hernández, Michael Tobler, Gary Voelker 2004;Agorreta et al., 2013).The mountains at PDM serve as a filter for primary and secondary Neartic and Neotropical fishes (Obregon-Barboza et al., 1994), frequently setting distribution limits and influencing phenotypic differences (Hulsey et al., 2010).For P. mexicana, unique haplotypes exist north of the PDM, which may be indicative of such isolation.In contrast, our results for P. sphenops and P. limantouri, which exhibited shared haplotypes on either side of the TMVB, suggested that the PDM is not a strong isolating barrier for these species at this time.However, the presence of P. sphenops outside of the described range (Miller et al., 2005), especially north of the TMVB, may also be a result of human introductions (Espinosa Perez, Ramírez, 2015).
The Balsas River system is the largest hydrological system on the Pacific slope of Mexico and houses many endemic species.The high rate of endemism in the Balsas Depression (BD) is a result of the isolation due to the formation of the Trans Mexican Volcanic Belt and Sierra Madre del Sur (Ferrusquía-Villafranca, 1998;Ferrari et al., 2000).The Balsas River system has been linked to events of dispersals and colonization by Neartic fish from the Mesa Central, a high plateau in central Mexico (Domínguez-Domínguez et al., 2010), a pattern coinciding in fish parasites (Martínez-Aquino et al., 2014), or with isolation of populations due the uplift of the TMVB (Pérez-Rodríguez et al., 2015).In our case, the Balsas River Valley is a potential area of transition for the northern expansion of P. nelsoni and P. sphenops.For P. nelsoni, the southern distribution of our sampling in Chiapas showed higher genetic diversity, potentially in relation to older populations, and for P. sphenops, the BD formed part of the historical documented range (Miller et al., 2005).This idea is supported by the shared populations or sister taxa in surrounding areas of the BD for Mollienesia and other species and coincides with the most recent geological isolation of this area by Plio-Pleistocene volcanism (Marshall, Liebherr, 2000).Future fine scale sampling may reveal additional patterns in the BD, as the Sierra de Taxco divides the Balsas basin into western and eastern units (Cabral-Cano et al., 2000).Interestingly, this ichthyogeographical province (Miller, Smith, 1986;Miller et al., 2005) houses the native endemic species P. maylandi of the subgenus Mollienesia, but sampling efforts in and around the type locality (Meyer, 1983) were unsuccessful, instead finding only P. nelsoni and P. sphenops.
Finally, the Isthmus of Tehuantepec is a low elevation area recognized as a strong faunal barrier that emerged during the Pliocene and Pleistocene (Campbell, 1999;Marshall, Liebherr, 2000;Mulcahy et al., 2006).This potential barrier does not appear to have played a significant role in shaping inter-specific variation in Mollienesia taxa (Fig. 2).However, our results for P. mexicana suggest that the Isthmus may be a moderate barrier, with one sub-clade being distributed almost entirely south of the Isthmus of Tehuantepec (and this sub-clade was sister to Guatemalan samples), and another being distributed entirely north of the Isthmus of Tehuantepec (Fig. 5).The same pattern did not hold for either P. sphenops or P. nelsoni, as both species had shared haplotypes occurring across the Isthmus (Figs. 3, 4).Our results for these species contrast with a recent study of Poeciliopsis pleurospilus and P. gracilis which, based on ddRADseq data, showed that the Isthmus of Tehuantepec was a clear geographic barrier between them (Ward et al., 2022).
In summary, the phylogenetic patterns of Mexican species in the subgenus Mollienesia demonstrate several independent invasions from Central America and subsequent diversification of species associated with isolation, both because of major Mexican molly phylogeography physiogeographic barriers and ecological selection.The phylogeographic patterns observed within Mexican species either indicated strong genetic structuring across a physiographic barrier (P.nelsoni and P. mexicana) or a lack of genetic structuring (P.sphenops, P. butleri, P. limantouri, and P. sulphuraria/thermalis), indicating that different species may vary significantly in their dispersal abilities.In fact, P. sulphuraria and P. thermalis are so closely related that they may represent the same species, but a formal taxonomic revision is required to determine whether P. sulphuraria should be treated as a synonym to P. thermalis (Palacios et al., 2013).Our analyses also indicated broader distributions than currently described for several species.This study also sheds insight into the taxonomic status of Pacific populations of P. sphenops, which was morphologically hypothesized as distinct based on the scale count around the peduncle (Miller et al., 2005).Since we find lack of genetic differentiation across the species' distribution, these morphological differences are likely indicative of geographic variation in phenotypes and not species differences.

FIGURE 1 |
FIGURE 1 | Sampling localities of species in the subgenus Mollienesia in Mexico and the main physiographic barriers throughout the country.The main physiographic barriers in Mexico from north to south are Northwestern Plains and Sierras, Sierra Madre Occidental, Sierra Madre Oriental, Gulf Coast Plain, Trans Mexican Volcanic Belt, Balsas Depression, Sierra Madre del Sur, Isthmus of Tehuantepec, and Sierra Madre de Chiapas.

FIGURE 2 |
FIGURE 2 | Bayesian tree from the MrBayes partitioned analysis of Poecilia spp.for two mitochondrial genes (Cyt b and ND2, 2187 base pairs) and one nuclear (RAG1, 1561 base pairs) rooted with other poeciliid outgroups.Species names in black pertain to species outside of Mexico and species names in red are species found within Mexico.Species labels represent slope: (A) = Atlantic, (P) = Pacific, and (B) = Bi-coastal.Nodal support shown (left to right; respectively): Bayesian Posterior Probabilities followed by RAxML bootstrap support values.Asterisks denote nodal support of 95% or above for the two methods, and a single asterisk at a node indicates support values of 95% or above for both methods.Nodes with no values present either had low values or were of little interest for this study.The capital letters at the end of each sample represents the state of origin in Mexico, from North to South: T = Tamaulipas, NVL = Nuevo León, V = Veracruz, H = Hidalgo, M = Michoacan, G = Guerrero, O = Oaxaca, Tb = Tabasco, C = Chiapas.
FIGURE 3 | Cytochrome b (1,140 bp) mitochondrial gene Bayesian phylogeny, parsimony haplotype network, and sampling distribution of Poecilia sphenops (aqua-Atlantic North of the Trans Mexican Volcanic Belt, pink-Balsas River Drainage, purple-Atlantic North of the Isthmus of Tehuantepec, red-Pacific North of the Isthmus of Tehuantepec, orange-Atlantic South of the Isthmus of Tehuantepec, and yellow-Pacific South of the Trans Mexican Volcanic Belt) across geographic barriers along the coasts of Mexico.The phylogeny has Bayesian posterior values followed by bootstrap values with asterisks representing support of 95% or above.The Parsimony network values correspond to the haplotype values and are colored according by geographic locations separated by barriers; circle sizes correspond to the number of individuals with that haplotype (larger circles reflect more individuals), and black circles indicate unsampled haplotypes.The capital letters at the end of each sample represents the state of origin in Mexico, from North to South: V = Veracruz, H = Hidalgo, M = Michoacan, G = Guerrero, O = Oaxaca, C = Chiapas.
FIGURE 4 | Cytochrome b (1,140 bp) mitochondrial gene Bayesian phylogeny, parsimony haplotype network, and sampling distribution of the Pacific sister taxa Poecilia butleri (blue-North of the Trans Mexican Volcanic Belt) and P. nelsoni (lime green-South of the Trans Mexican Volcanic Belt, pink-Balsas River Drainage, red-North of the Isthmus of Tehuantepec, and yellow-South of the Trans Mexican Volcanic Belt) across geographic barriers along the coast of Mexico.The phylogeny has Bayesian posterior values followed by bootstrap values with asterisks representing support of 95% or above.The Parsimony network values correspond to the haplotype values and are colored according by geographic locations separated by barriers; circle sizes correspond to the number of individuals with that haplotype (larger circles reflect more individuals), and black cirles indicate unsampled haplotypes.The capital letters at the end of each sample represents the state of origin in Mexico, from North to South: S = Sinaloa, N = Nayarit, J = Jalisco, Cl = Colima, M = Michoacan, G = Guerrero, O = Oaxaca, C = Chiapas.

FIGURE 5 |
FIGURE 5 | Cytochrome b (1,140 bp) mitochondrial gene Bayesian phylogeny, parsimony haplotype network, and sampling distribution of Atlantic taxa Poecilia limantouri (Turqouise-North of the Trans Mexican Volcanic Belt, forest green-South of the Trans Mexican Belt), P. sulphuraria/P.thermalis (yellow-South of the Isthmus of Tehuantepec), and P. mexicana (baby blue-North of the Trans Mexican Volcanic Belt, light green-South of the Trans Mexican Volcanic Belt, purple-North of the Isthmus of Tehuantepec, orange-South of the Isthmus of Tehuantepec) across geographic barriers along the Atlantic coast of Mexico.The phylogeny has Bayesian posterior values followed by bootstrap values with asterisks representing support of 95% or above.The Parsimony network values correspond to the haplotype values and are colored according by geographic locations separated by barriers; circle sizes correspond to the number of individuals with that haplotype (larger circles reflect more individuals), and black circles indicate unsampled haplotypes.The capital letters at the end of each sample represents the state of origin in Mexico, from North to South: NVL = Nuevo Leon, T = Tamaulipas, V = Veracruz, H = Hidalgo, Tb = Tabasco, C = Chiapas.

TABLE 1 |
List of species from the Poecilia sphenops and P. mexicana species complexes present in Mexico.The table provides information about the slope on which a species occurs, a broad description of the distribution, the endemism status in Mexico, and phylogeographic breaks uncovered in this study.Asterisks indicate the inclusion of a species in phylogenetic analyses; crosses indicate focal species used for phylogeographic analyses based on Cyt b.