Phylogeography of the Mesa Silverside fish Chirostoma jordani (Woolman, 1894) throughout the Mexican Plateau

Sample collection and molecular data

We obtained 170 specimens of C. jordani from 33 localities across the species’ distributional range (Table 1, Fig. 1), including the only two known northern populations: one from the Bolaños River and the other from the upper Mezquital. We captured fishes with the permission of local authorities, using electrofishing and seine netting methods, and euthanized them with tricaine-mesylate (MS-222). In addition, some specimens were obtained by local fishermen. The care and use of the animals complied with animal welfare laws, guidelines and policies, as approved by SEMARNAT-SGA/DGVS/2009/19, SEMACCDET-OS-0084/2019 and PPF/DGOPA-014/20. We preserved fin clips tissue taken from each specimen in 96% ethanol. The whole fish were fixed in 5% formalin, preserved in 70% ethanol, and identified following Barbour (1973b) and Miller, Minckley & Norris (2005). Tissue samples and specimens were deposited in the Ichthyological Collection of the Universidad Michoacana de San Nicolas de Hidalgo (CPUM).

We isolated genomic DNA using the conventional proteinase K/phenol/chloroform protocol (Hillis, Moritz & Mable, 1996). Polymerase chain reaction (PCR) was used to amplify the mitochondrial Cytochrome b locus (Cytb: 1,140 bp) using the primers Glud-G (Palumbi, 1996) and H16460 (Perdices et al., 2002). A fragment of the hypervariable control region (D-loop: 350 bp) also was amplified using the primers RCA and RCE (Anderson et al., 1981; Kim et al., 1999) and, a fragment of the first intron of the ribosomal protein S7 (S7: 629 bp) also was amplified using the primers S7RPEX1F and S7RPEX1R (Chow & Hazama, 1998). The PCR cycling conditions and primers sequences are described in Table S1. The PCR products were visualized on a 1.5% agarose gel and amplicons were submitted to MACROGEN Korea for sequencing.

The chromatograms of the recovered sequences were manually examined to eliminate potential sequencing errors. We performed the alignments manually using the software Geneious v.4.8.5. We translated the Cytb alignments to amino acids to verify the absence of stop codons along the sequence. We resolved the heterozygous sites for the S7 sequences using the algorithm provided by PHASE 2.0 (Stephens, Smith & Donnelly, 2001), as implemented in DnaSP v. 5.10 (Librado & Rozas, 2009), employing the default parameters. The non-recombination test of S7 was conducted using the software DnaSP v. 5.10 (Librado & Rozas, 2009). Additionally, we also gathered sequences data from several related species including Chirostoma estor, Chirostoma humboldtianum, Chirostoma chapalae, Chirostoma sphyraena, and Chirostoma attenuatum. These species were used as outgroups in the gene tree analyses. The DNA sequences were deposited in NCBI databases (https://www.ncbi.nlm.nih.gov/), accession numbers are included in Table S2. For the phylogenetic analyses, the best-fit model for each marker was obtained in PartitionFinder v.1.1.0 (Lanfear et al., 2012) using the Akaike Information Criterion (AIC).

Gene trees, genetic distance and haplotype networks

We inferred gene trees of haplotype data files of both mtDNA (Cytb and D-loop) and nDNA (S7) using both, maximum likelihood (ML) and Bayesian Inference (BI) approaches. The ML was run in RAxML (Stamatakis, Hoover & Rougemont, 2008; Silvestro & Michalak, 2012; Stamatakis, 2014) with the AICc-select model with three gene partitions. Bayesian analyses were run in MrBayes 3.2.2 (Ronquist et al., 2012). Two independent runs were implemented with four MCMCs and 10,000,000 generations, sampling every 100 trees. Chain convergence was verified by a suitable effective sample size (ESSs >200) for all parameters in Tracer package v. 1.7, discarding 10% of generations as burn-in (Rambaut et al., 2018; Sahlin, 2011). Additionally, we estimated the uncorrelated p genetic distance based on both mtDNA and nDNA using “APE” v.5.8 (Paradis, Claude & Strimmer, 2004) package implemented in R.

To examine the geographic distribution of haplotypes, we constructed a haplotype network for each independent locus in HaploView v. 4.2 (Barrett et al., 2005), based on an unrooted ML tree obtained in RAxML, through the CIPRES Science Gateway v.3.3 (Miller et al., 2015).

Genetic structure

To determine the amount of genetic variation partitioned within and among populations, an analysis of molecular variance (AMOVA) was performed in ARLEQUIN (Excoffier, Smouse & Quattro, 1992; Excoffier & Lischer, 2010) for the three separate loci. This analysis was run with significance levels set at α = 0.05 and 10,000 random permutations.

We used Hierarchical Bayesian Analysis of Populations Structure (hBAPS) (Cheng et al., 2013), to examine the genetic population structure, using the Single Nucleotide Polymorphism (SNPs) matrix from the mtDNA and nDNA sequences. The analysis was carried out in the RhierBAPS (rBAPS) package implemented in R language (Cheng et al., 2013; Hill et al., 2019). The initial number of clusters assigned was equal to the number of populations with a length of run = 1,000,000. To plot the groups obtained, we used a guide tree result estimated using a pml: Likelihood in package “phangorn” (Schliep, 2011) implemented in R. The bar-plot graphics of probability assignment of rBAPS results were generated in R package. Additionally, we calculated the Φct, Φst and Φsc values in ARLEQUIN (Excoffier, Smouse & Quattro, 1992; Excoffier & Lischer, 2010) for the main clusters recovered in hBASP to mtDNA and nDNA.

Time calibration

To estimate the divergence time within C. jordani we concatenated mtDNA and nDNA haplotypes into a single data set. The time calibration analysis was performed using BEAST 2 package (Bouckaert et al., 2014), using a coalescent prior (Bayesian skyline plot). We included samples of related species as described in the gene tree analysis. We used an uncorrelated lognormal clock. Due to a limited fossil record for the group, the molecular clock was calibrated with Cytb rate of 0.005 substitutions/site/Mya estimated for Atheriniformes fishes (Campanella et al., 2015; Unmack, Allen & Johnson, 2013). We estimated the mutation rate of D-loop and S7 relative to Cytb, using a normal prior (0.005 ± 0.001). The best-fit models of nucleotide evolution were used as estimated in PartitionFinder. The assays were performed in BEAST 2 implemented on the web server CIPRES Science Gateway v. 3.3 (Miller et al., 2015). The analysis was run for 100,000,000 generations sampling each 1,000 generations. Ten percent of the runs were discarded as burn-in. Two independent runs were performed and both files were combined using LogCombiner v.2 (Bouckaert et al., 2014). Posterior probability density of the combined tree file was summarized with TreeAnnotator v. 2 (Helfrich et al., 2018).

Molecular data

One hundred seventy specimens were collected from 33 localities from 14 biogeographic regions (sensu Domínguez-Domínguez et al., 2010) in the Mexican Plateau (Fig. 1; Table 1). One hundred forty-five sequences were obtained for the mitochondrial Cytb locus of 1,000 bp of length. One hundred seventy sequences were obtained for the mitochondrial D-loop locus with 333 bp. For the S7 nuclear locus, 99 sequences were obtained with 629 bp, and all heterozygous sites found in the S7 sequences were successfully resolved for SNP variation (phase threshold value >85%), producing 146 genotypes (Table 1). The best-fit models of nucleotide evolution were as follows, for Cytb was TrN + G, GTR +I +G for D-loop and GTR+G for the S7.

Gene trees, genetic distance and haplotype networks

The mtDNA gene tree using BI recovered two clades (Fig. 2). Clade I (Fig. 1) is not supported by BI but is supported by ML and includes haplotypes from La Vega (Ameca basin), Belem del Refugio, Nochistlan, Lagos de Moreno, and Tepatitlan in the Verde River (Verde-Santiago Basin), as well as the Guadalupe Victoria dam (Mezquital basin). Clade II (Fig. 1) is well supported by BI (posterior probabilities (pp) = 1) but not by ML, and includes haplotypes from Cazones, Chapala, Cotija, Cuitzeo, Bolaños-Santiago, Magdalena, Lerma, Pánuco, Valle de México, El Chichimeco, La Paz, Ojuelos, Río el Chilerillo, and San Nicolas de las Flores (Verde-Santiago basins) (Fig. 2). The ML estimation also recovered two clades, with the same population distribution as in the BI analysis. Clade I is supported with a bootstrap of 87, but Clade II is not supported by ML (Fig. 2). The S7 haplotype data gene tree recovered a basal polytomy among populations without geographical clustering in both BI and ML approaches (Fig. S1). The uncorrected P-distance between Clade I and Clade II was estimated at 3% using mtDNA sequences (Fig. S2) and 1% with nDNA (Fig. S3).

The haplotype network analysis of the mitochondrial locus reveals the presence of two distinct haplogroups. Within Haplogroup 1, internal differentiation is observed, with four geographically distinct groups: two distributed in different tributaries of the Verde-Santiago drainage, separated by four to seven mutation steps from the nearest haplotype; one in the Ameca drainage, separated by three to five mutation steps; and one in the Mezquital drainage, separated by ten to eleven mutation steps (Figs. 1 and 2A–2B). Haplogroup 2 clusters the remaining samples from western, central, and eastern drainages, including the Verde River (within the Santiago drainage), as well as samples ranging from Magdalena in the west (number 13 in Fig. 1) to Cazones in the east (number 3 in Fig. 1). The S7 network does not show distinct haplogroup formation, although the Mezquital and Ameca samples display genetic distinctiveness, with some mixing observed between the Verde and Mezquital populations. The remaining locations share common haplotypes (Fig. 3C).

Genetic structure

The AMOVA analysis of the mtDNA, without a priori groupings, showed a high and significative Φst value (Φst = 0.51), rejecting the null hypothesis of panmixia. The majority of the total variation (50.61%) was explained by differences between populations (Φst>Φsc). In contrast, for the nDNA sequences, the majority of the variation (67.1%) was explained within populations, with only 35.52% explained between populations (Φst<Φsc). The Φst value differed significantly from 0 (Φst = 0.329), again rejecting the null hypothesis (Table 2).

The hBAPS analysis for mtDNA partitioned the genetic variation in six clusters. Genetic cluster (K1) included samples from Cuitzeo, Bolaños, Cotija, Magdalena, Middle Lerma, Chapala, Santiago, Verde Santiago, Panuco and Valle de México. The second genetic cluster (K2) grouped samples from Cazones, Lower Lerma, and Verde Santiago. The third genetic cluster (K3) grouped samples from the Ameca and, Verde Santiago. The fourth (K4) and fifth (K5) genetic clusters grouped samples from Verde Santiago. Finally, the sixth genetic cluster (K6) grouped samples from Mezquital (Figs. 2 and 4). The F-statistics for this arrangement were: Φct = 0.36, Φst = 0.55 and Φsc = 0.29 (Φct>Φsc). For nDNA, we recovered three genetic clusters, but with extensive admixture between geographical locations. The first cluster (K1) grouped samples from Chapala, Cuitzeo, Middle and Low Lerma, Magdalena and Verde Santiago Basin. The second genetic cluster (K2) grouped samples from Ameca, Bolaños-Santiago, Cazones, Chapala, Cuitzeo, Cotija, Lower Lerma, Magdalena, Middle Lerma, Verde Santiago and Panuco basin, and the third genetic cluster (K3) grouped samples of Mezquital and Verde Santiago (Nochistlan) basin (Fig. 5). The F-statistics for this arrangement were: Φct = 0.17, Φst = 0.38 and Φsc = 0.25.

Time calibration analysis

The divergence time analysis showed the split between clade I and clade II dated in the Pleistocene ca. 1.4 Mya (95% HDP 0.87−2.1 Mya) (Fig. 6).

Split of the two main clades

Our results support the existence of two main clades within C. jordani (Fig. 2). The split was dated ca 1.4 Mya (Fig. 6). Clade I is confined to western (Ameca, Verde-Santiago) and northern basins (Mezquital), while clade II is widely distributed in western, central, and eastern basins on the Mexican Plateau (Figs. 1 and 2), showing mixtures of haplotypes between regions, except for the Valle de Mexico, Panuco and Cazones samples for D-loop (Fig. 3A). The biogeographic scenario that explains the cladogenetic events in C. jordani is complex, mainly because the mixture of populations belonging to the two haplogroups in the Verde River. This area has been geologically active in the last 3 Mya. We hypothesize that the isolation between Clade I and Clade II and the mixture of populations belonging to different haplogroups, has been promoted by the high tectonic and volcanic activity in the “triple junction area” between 3 and 1 Mya, which allowed for the isolation and reconnection of the upper Ameca, Verde, Lerma Rivers and Chapala Lake regions (Rosas-Elguera & Urrutia-Fucugauchi, 1998; Nieto-Samaniego et al., 1999). The influence of the geological activity of the “triple junction” in the evolution patterns of the Mexican ichthyofauna has been extensively discussed by other authors (Barbour, 1973b; Smith, Cavender & Miller, 1975; Miller, Minckley & Norris, 2005; Domínguez-Domínguez et al., 2010; Pérez-Rodríguez et al., 2016; Beltrán-López et al., 2018; Beltrán-López et al., 2021). The climatic instability during the early Pleistocene, a period characterized by intercalated humid and dry periods, has been thought to have had an important influence on the evolutionary and demographic history of the Mexican ichthyofauna (Barbour, 1973b; Mateos, Sanjur & Vrijenhoek, 2002; Domínguez-Domínguez et al., 2010; Pérez-Rodríguez et al., 2009a; Pérez-Rodríguez et al., 2009b; Pérez-Rodríguez et al., 2015; Beltrán-López et al., 2018; Betancourt-Resendes, Pérez-Rodríguez & Domínguez-Domínguez, 2018; Betancourt-Resendes et al., 2019; García-Martínez et al., 2015; García-Martínez et al., 2020). The climatic scenario may have influenced the interchange of fauna in the tributaries of the Verde and Lerma Rivers, promoting the connection of headwater tributaries in wet periods, and the mixture of the Verde population in both haplogroups,

Phylogeographic patterns of clade I

Our phylogeographic and population genetic results support well-structured genetic groups with geographic correspondence within Clade I/Haplogroup 1 based on mtDNA and partially for nDNA (Figs. 2, 3, 4 and 5). The TRMCA was dated ca. 1 Mya (Fig. 6). Within this clade, one group occurred in the north (Mezquital Basin in Mesa del Norte) and three in the west (Ameca and Verde-Santiago basins in central Mexico).

The biogeographic relationship between the Mezquital Basin and the central basins of Mexico has been widely recognized. The genus Characodon (Goodeidae) and the species Moxostoma milleri (Catostomidae), are distributed in the Mezquital River basin and have their closest relatives distributed in central Mexico (Doadrio & Domínguez, 2004; Domínguez-Domínguez et al., 2010; Pérez-Rodríguez et al., 2015). However, the phylogenetic relationships between these groups show older divergence times than Chirostoma (<1Mya), being ca. 4.5 Mya for the split of M. milleri, and ca. 15.5 Mya for the split of Characodon. Although is not clear how the Mezquital River was connected to the central Mexican basins in recent geological times, according to the peripheral position of the Mezquital samples in the three haplotype networks, a plausible isolation scenario could be the range extension of the MRCA from the Verde-Santiago drainage to the Mezquital River, with a posterior vicariant event that generated the allopatric pattern. The likely connection between the Mezquital River and central Mexico could be through the Guadiana Valley, a high plain in the upper Mezquital that border in its southern portion with the rivers that drain to the Santiago River. The accumulation of Pleistocene alluvial deposits within the Guadiana Valley, in synergy with the volcanic activity (Albritton, 1958; Aranda-Gómez et al., 2018) and the climatic changes that occurred during Pleistocene, could have caused the isolation of the MRCA of C. jordani from Mezquital drainage.

The taxonomic status of C. jordani from the Mezquital drainage has been highly discussed. In morphological analyses, the Mezquital samples show extensive morphological polymorphism when compared with specimens of C. jordani collected from throughout the entire distribution of the species, and has therefore, largely been considered a synonym of C. jordani (Barbour, 1973a). Even in the Integrated Taxonomic Information System-ITIS (Last accessed 10/24/2023, http://www.itis.gov) and FishBase (Last accessed 10/24/2023, https://www.fishbase.se/summary/Chirostoma-jordani.html), it appears as a synonym of C. jordani. However, Miller, Minckley & Norris (2005) recognized C. mezquital as a valid species based on characters associated with body height measurements. The results presented herein support the genetic distinctiveness of the Mezquital samples, forming an independent group in haplotype network, and hBAPS analyses (Figs. 2–5). However, other populations, including the Ameca and Verde River populations, also show genetic differences. According to the results presented in this study, we propose an exhaustive integrative taxonomic study to elucidate the taxonomic status of C. jordani population from Mezquital and other genetically differentiated populations.

The samples collected in the Verde River basin show three genetically segregated groups inhabiting this basin. Two groups were clustered into Clade I and the other group of samples was clustered in Clade II (Figs. 2–5). As previous studies mentioned, the existence of genetically differentiated populations within the Verde River could be the result of the existence of temporally and spatially independent events of colonization in the Verde River (Figs. 2–5), as was previously proposed for Poeciliopsis infans, Goodea atripinnis, Algansea tincella and Notropis calientis (Domínguez-Domínguez et al., 2008a; Domínguez-Domínguez et al., 2009; Pérez-Rodríguez et al., 2009b; Beltrán-López et al., 2018; Beltrán-López et al., 2021). Hence, we hypothesize that the first colonization event to the Verde River was through the postulated connection of the Verde Paleoriver and Chapala Paleolake (Smith, Cavender & Miller, 1975; Barbour, 1973a; Miller & Smith, 1986; Aranda-Gómez, Henry & Luhr, 2000; Miller, Minckley & Norris, 2005). A second colonization event seems to have occurred through a river capture event between the Middle Lerma and Verde River. The location of Lagos de Moreno, Nochistlan and Belen del Refugio is geographically close to the headwaters of the Turbio River, a tributary of the Lerma River. In this case, the intermontane plains that separate both river systems could function as corridors during wet and flood periods. This hypothetical connection has been recognized previously as a dispersal route for the ancestor of the species complex belonging to the atherinopsid silversides and other species (Barbour, 1973b; Bloom et al., 2012). A third invasion of the Verde River seems to be also related to stream captures of the Lerma and Verde Rivers; these events appear to have been recent enough to prevent haplotype sorting in any of the three sampled loci (Figs. 2–5).

The isolation and cladogenesis of the samples of Chirostoma from the Ameca basin are not surprising, since seven endemic freshwater fishes that occurred in this river basin have their closest relatives in Central Mexico drainages, like the leuciscids Algansea ameca, Notropis amecae and Yuriria amatlana, the goodeids Zoogoneticus tequila, Allodontichthys polylepis, Xenotoca doadrioi, Skiffia francesae and Allotoca goslinei, and the catostomid Moxostoma mascotae (Miller & Smith, 1986; López-López & Paulo-Maya, 2001; Doadrio & Domínguez, 2004; Domínguez-Domínguez, Pompa-Domínguez & Doadrio, 2007a; Domínguez-Domínguez, Pompa-Domínguez & Doadrio, 2007b; Domínguez-Domínguez et al., 2008a; Domínguez-Domínguez et al., 2009; Domínguez-Domínguez et al., 2010; Dominguez-Dominguez, Doadrio & Perez-Ponce de Leon, 2006; Piller et al., 2015; Pérez-Rodríguez et al., 2009a; Pérez-Rodríguez et al., 2009b; Pérez-Rodríguez et al., 2016; Beltrán-López et al., 2018). The high endemicity of the Ameca River seems to be related to a long history of volcanism and tectonism at different geological times, generating the isolation, connection and compartmentalization of the Ameca River over time, which includes several surrounding drainages (Rosas-Elguera & Urrutia-Fucugauchi, 1998; Garduño & Tibaldi, 1991). Our results support a recent isolation event of the Ameca River, since we date the cladogenesis of the Ameca group at ca. 0.90 Mya (Fig. 6). Two recent hydrographic connections and interchange of fauna between Ameca River and the contiguous basins have been postulated with the first through the Atotonilco, San Marcos and Zacoalco-Ameca Paleolakes area (Smith, Cavender & Miller, 1975), with a subsequent isolation event, occurred less than ca 1 Mya. This isolation seems to be promoted by Pleistocene volcanism and the intense tectonic activity of the so-called triple junction (Rosas-Elguera & Urrutia-Fucugauchi, 1998; Garduño & Tibaldi, 1991), as has also been postulated for the isolation of P. infans, Ameca splendens, Xenotoca melanosoma and Z. purhepechus populations in the Ameca River (Domínguez-Domínguez et al., 2008a; Beltrán-López et al., 2018). A second connection event was proposed for P. infans between the Ameca River and the Verde-Santiago River drainage, suggesting that these basins were connected until very recent geological time through stream capture of the Ameca and Verde Rivers, which was facilitated by the volcanism in the Tepic-Zacoalco graben (Mateos, Sanjur & Vrijenhoek, 2002; Beltrán-López et al., 2018). However, we cannot support either of these two scenarios, due to the lack of samples of Chirostoma from the Sayula-San Marcos-Atotonilco lakes area, despite the high sample effort conducted.

Phylogeographic patterns of Clade II

The Clade II, shows a widely distributed group in the central and eastern drainages of central Mexico (in an area > 130,000 km²), from Bolaños River in the Santiago Drainage in the Northwest, Ameca and Magdalena in the Southwest, to Cazones to the East, with a high haplotype mixture (Figs. 3 and 5). The group has been evolving in central Mexico over the last ∼1 Myr. This genetic cohesiveness of widely distributed species in central Mexico is not common, since most of the co-distributed species already studied are structured in at least one of the locations where haplotypes of the C. jordani central group are distributed, as is the case of Zoogoneticus quitzeoensis (Domínguez-Domínguez et al., 2008a), Notropis calientis (Domínguez-Domínguez et al., 2009), Moxostoma austrinum (Pérez-Rodríguez et al., 2016), Algansea spp (Pérez-Rodríguez et al., 2009a; Pérez-Rodríguez et al., 2009b), and Poeciliopsis infans (Mateos, Sanjur & Vrijenhoek, 2002; Beltrán-López et al., 2018).

The genetic homogeneity could be related to ancient connectivity between the central basins across several periods during the Pleistocene, that have been widely discussed in goodeids (Doadrio & Domínguez, 2004; Domínguez-Domínguez et al., 2010; Beltrán-López et al., 2021), leuciscids (Schönhuth, 2002; Domínguez-Domínguez, Pompa-Domínguez & Doadrio, 2007a; Domínguez-Domínguez, Pompa-Domínguez & Doadrio, 2007b; Pérez-Rodríguez et al., 2009a; Pérez-Rodríguez et al., 2009b; García-Andrade et al., 2021) and poecilids (Mateos, Sanjur & Vrijenhoek, 2002; Beltrán-López et al., 2018), as well as to the dispersal capacity of the species. In fishes, this capacity is highly dependent on body shape (Sfakiotakis, Lane & Davies, 1999; Triantafyllou, Triantafyllou & Yue, 2000; Langerhans & Reznick, 2009), whereas the survival rate depends on plasticity and adaptability for establishment in a new habitat (Seehausen & Wagner, 2014). In particular, C. jordani exhibits a high degree of morphological variation (Barbour, 1973b), part of which is associated with divergent habitats, and differential selective pressures (Foster, Bower & Piller, 2015). Moreover, its small size and ecological versatility (Barbour, 1973a) gives C. jordani the capacity to widely disperse throughout the lakes, rivers, and streams of central Mexico, and also to survive in newly invaded environments, taking advantage of possible historical dispersal routes in currently isolated basins (Domínguez-Domínguez et al., 2008a; Pérez-Rodríguez et al., 2009a; Pérez-Rodríguez et al., 2009b; Beltrán-López et al., 2018).