Background {#Sec1}
==========

The structure of current North American biodiversity is related to geological (tectonic and volcanic) and climatic changes that occurred mainly during the Neogene (*ca.*  33 Million years ago (Mya) and Quaternary (*ca.*  past 2.5 Mya) \[[@CR1], [@CR2]\].  The Sierra Madre Occidental (SMOC) mountain range in Mexico is an area of high faunal and floristic endemism, and is considered an important biogeographic corridor \[[@CR3]\].  The formation of the SMOC during the Oligocene (*ca.*  33-23 Mya) until the present \[[@CR2], [@CR4]\] is an important factor influencing the evolutionary history of several taxa, particularly organisms with low vagility such as freshwater fish species \[[@CR3], [@CR5]--[@CR12]\]. 

Although the biodiversity of northwestern Mexico is mainly related to tecto-volcanic activity of the Tertiary, the SMOC is also considered a Pleistocenic refuge, and associated with the expansion and contraction of the distribution of numerous species \[[@CR13]\].  The effects of climate changes on epicontinental waters has significantly affected the distribution of the fish fauna of this region \[[@CR3]\].  The evolutionary history of the sucker fish family Catostomidae in Mexico may have been influenced by these processes, especially in the case of the Catostominae subfamily, the predominant fish group of the family in Mexico \[[@CR14]\]. 

Within the Catostominae, two closely related species of the genus *Pantosteus* Cope & Yarrow, 1875; inhabit northern Mexico including the SMOC: *Pantosteus plebeius* (Baird & Girard, 1854) and *Pantosteus nebuliferus* (Garman, 1881) \[[@CR14], [@CR15]\]. * Pantosteus nebuliferus* has a restricted range, and is endemic to the endorheic Nazas and Aguanaval drainage basins.  In contrast, *Pantosteus plebeius* is widespread and occurs in several drainage basins across the SMOC: the Mezquital, Piaxtla, Fuerte, and Yaqui basins; the Central Guzman hydrographic system (Santa Maria, Casas Grandes, and Del Carmen basins); the Upper Rio Grande basin (Conchos River in Mexico and Rio Grande River in USA); and the Mimbres basin in New Mexico, USA (Fig.  [1](#Fig1){ref-type="fig"}) \[[@CR16]\]. * Pantosteus plebeius* is also distributed in the Gila River of the Colorado basin, but the origin of this population is uncertain, some consider it to have been artificially introduced \[[@CR17], [@CR18]\], while others suggest this population is the result of stream capture from the Mimbres River \[[@CR19]\].Fig.  1Distribution of *Pantosteus plebeius-nebuliferus* species-group and sampling localities.  The colours correspond to the eight main basins and numbers to the localities: (1) Santa Gertrudis, (2) Arroyo Las Bayas, (3) Pino Suarez, (4) La Barranca, (5) Puente Mimbres, (6) El Cuarto, (7) Peñón Blanco, (8) Covadonga, (9) Arroyo Torreones, (10) El Peñasco, (11) El Olote, (12) Atotonilco, (13) Las Vegas, (14) Fuerte, (15) Oteros, (16) Bocoyna, (17) Conchos, (18) Ureyna, (19) South Fork Palomas Creek, (20) Escalariado, (21) Casas Grandes, (22) Ignacio Zaragoza, (23) Buenaventura and (24) Santa Clara

As with other *Pantosteus* species in Mexico, *P. plebeius* is believed to have colonized several drainage basins during the postglacial period \[[@CR20]\].  Dispersal during anastomosis of rivers, followed by vicariance after river isolation, has been hypothesized as the main factor influencing genetic divergence among *Pantosteus* populations.  It has been suggested that the populations of some drainage basins could represent one or more independent evolutionary lineages, such as those of the Mezquital, Piaxtla, Fuerte, and upper Rio Grande basin, each displaying significant genetic structure compared with Mimbres and Gila River populations \[[@CR21], [@CR22]\].  Given the complex evolutionary history of *P. plebeius* it has been suggested that this lineage may represent several distinct species, one or more of which may be synonymous with *P. nebuliferus* \[[@CR23]\]. 

*Campostoma ornatum* (Girard, 1856) and *Codoma ornata (*Girard, 1856) are codistributed with the *Pantosteus plebeius-nebuliferus*c species-group throughout its range.  Phylogeographic studies of these species have shown repeated isolation scenarios followed by merge or convergence in response to range expansion (reticulated biogeographic history) \[[@CR3], [@CR11]\].  This pattern was related to repeated events of dispersal and isolation associated with tecto-volcanic activities occurring in the SMOC since the early Pliocene, as well as to climatic fluctuations during glacial and interglacial periods.  Population structure in *Campostoma ornatum* and *Codoma ornata* is reflected in high genetic differentiations among populations inhabiting the drainage basin in the northern SMOC, including the Guzman system, and the Río Grande, Yaqui, Mayo, Fuerte, and Nazas basins \[[@CR3], [@CR11]\]. 

The aim of this study was to 1) examine the phylogeography and reconstruct the evolutionary history of the *Pantosteus plebeius-nebuliferus* species-group in Mexico; 2) to explore the genetic diversity and structure of the species-group; and 3) determine population divergence times using nuclear and mitochondrial markers.  We hypothesize that the genetic structure and distribution of the main clades of *Pantosteus plebeius-nebuliferus* are associated with the geological history of the SMOC and central Northern Mexico. 

Methods {#Sec2}
=======

Sampling and dna isolation {#Sec3}
--------------------------

We collected 87specimens with permission of relevant authorities in July 2013 from 20 localities in 8 drainage basins (Fig.  [1](#Fig1){ref-type="fig"} and Table [1](#Tab1){ref-type="table"}) using electrofishing and hand nets.  We included 6 samples from 4 localities from GenBank (GU937833, KJ441235 and KJ441237-KJ441240) for a total of 24 localities in 9 basins (localities 14, 18, 19 and 20 in Fig.  [1](#Fig1){ref-type="fig"}).Table 1Sampled populations of the *Pantosteus plebeius-nebuliferus* species-groupSpeciesBasin*cytb* individuals*GHI* individualsGenBank Accession*Pantosteus nebuliferus*Nazas2929*cytb*: MG203677-MG203705\
*GHI* allele 1: MG203712-MG203736, MG203784-MG203787\
*GHI* allele 2: MG203794-MG203820, MG203868-MG203869*Pantosteus plebeius*Mezquital128*cytb:* MG203619-MG203626, MG203669-MG203672\
*GHI* allele 1: MG203706-MG203711, MG203770-MF203771\
*GHI* allele 2: MG203788-MG203793, MG203860-MG203861Piaxtla44*cytb*: MG203673-MG203676\
*GHI allele 1*: MG203772-MG203775\
*GHI* allele 2: MG203862-MG203865Fuerte119*cytb*: MG203627-MG203636, KJ441240 \[[@CR61]\]\
*GHI* allele 1: MG203737-MG203744, MG203783\
*GHI* allele 2: MG203821-MG203829Uppers North R. Grande (Palomas)11*cytb*: KJ441237 \[[@CR61]\]\
*GHI*: GU937833 \[[@CR105]\]Upper South R. Grande (Conchos)22*cytb*: MG203637, KJ441239 \[[@CR58]\]\
*GHI* allele 1: MG203745-MG203746\
*GHI* allele 2: MG203830-MG203831Del Carmen109*cytb*: MG203658-MG203657, KJ441238 \[[@CR61]\]\
*GHI* allele 1: MG203760-MG203767, MG203776\
*GHI* allele 2: MG203849-MG203857Santa Maria88*cytb*: MG203650-MG203666,\
*GHI* allele 1: MG203756-MG203759, MG203777-MG203779, MG203782\
*GHI* allele 2: MG203842-MG203848, MG203867Casas Grandes1513*cytb*: MG203638-MG203649, MG203667-MG203668, KJ441235 \[[@CR61]\]\
*GHI* allele 1: MG203747-MG203755, MG203768-MG203769, MG203780-MG203781\
*GHI* allele 2: MG203832-MG203841, MG203858-MG203859, MG203866

Despite a considerable sampling effort, we were unable to obtain any specimens from the Aguanaval and Yaqui basins.  We obtained fin clips and preserved them preserved in absolute ethanol and stored them at -75 °C.  We fixed several specimens in formalin for identification and deposited in the Fish Collection of the Universidad Michoacana de San Nicolas de Hidalgo CPUM, Morelia, Mexico.  We released the remaining captured fish alive. 

We performed isolation of genomic DNA with BioSprint DNA Blood Kit QIAGEN according to the manufacturer's instructions.  We amplified the complete cytochrome *b* mitochondrial gene (*cytb*) using the primers GLuDG \[[@CR24]\] and H16460 \[[@CR25]\] and the nuclear 3^rd^ intron of the growth hormone copy I (*GHI*), using the primers GHI3F and GHI3R \[[@CR26]\]. 

We amplified DNA samples by Polymerase Chain Reaction (PCR), using 12 μl volume reactions with final concentrations of 0.2 μM of each primer, 0.25 mM of each dNTP, 1.5 M of MgCL~2~, and 1 U of *Taq* DNA Polymerase.  The PCR procedure consisted of 2 min at 95 °C followed by 35 cycles of 45 s at 94 °C for denaturation, 1.5 min at 46 °C and 52.2 ° C for primer alignment for *cytb* and *GHI* respectively, 2 min at 72 °C for synthesis, and a final extension of 5 min at 72 °C.  We quantified the PCR products by electrophoresis on 1.5% agarose gel and submitted to Macrogen Inc. (Korea) and htSEQ Inc (High-Throughput Sequencing, University of Washington, USA) for sequencing. 

Genetic variation {#Sec4}
-----------------

We estimated the genetic diversity within drainage basins for the two gene markers based on nucleotide diversity (π), haplotype diversity (h), and the proportion of segregating sites (Θ~S~), using the software DNAsp 5.0 \[[@CR27]\].  We obtained the genetic *p*-distances (D~P~) between populations using Mega 5.2 \[[@CR28]\].  A bootstrapping process was implemented with 1000 repetitions. 

We estimated the genetic structure at the interspecific and intraspecific level of the *Pantosteus plebeius-nebuliferus* species-group with fixation indices (Φ~ST~) for both molecular markers.  We used spatial analysis of molecular variance of geographically homogeneous K groups (SAMOVA) with SAMOVA 2.0 \[[@CR29]\] and AMOVA (for K=1 to identify genetic structure without *a priori* information and the genetic variance explained within population) with Arlequin 3.5.1.3 \[[@CR30]\] to define groups of populations (basins) that were maximally differentiated without constraints of their geographic distribution.  For both SAMOVA and AMOVA analysis we used 10 000 iterations from each of 100 random initial conditions.  We tested for each K value from 2 to 8 in SAMOVA.  We applied a Bonferroni correction \[[@CR31]\] to each *p*-value obtained in the paired test of genetic differentiation. 

Phylogenetic analysis and haplotype network reconstruction {#Sec5}
----------------------------------------------------------

We manually aligned DNA sequences in Mega 5.2 and examined them using chromatograms.  We phased *GHI* sequences with point mutations using DNAsp and applied a test of recombination using a coalescent-based Bayesian method (10 000 replicates) in the same software. 

We used two algorithms for phylogenetic reconstructions, using *Catostomus catostomus* as outgroup (GenBank Accession: AF454871 for *cytb* and GU937824 for *GHI*).  We conducted independent Maximum Likelihood (ML) analyses with RAxMLGUI 1.3.1 \[[@CR32], [@CR33]\] for both genes, performing 10 000 bootstrap repetitions and using the evolutionary substitution model estimated with PartitionFinder \[[@CR34]\].  These models were General Time-Reversible \[[@CR35]\] + gamma (GTR+G) for the complete *cytb* and Tamura-Nei \[[@CR36]\] + gamma (TrN+G) for *GHI* (see Additional file [1](#MOESM1){ref-type="media"}).  We conducted Bayesian phylogenetic (BI) reconstruction with the software MrBayes 3.2.6 \[[@CR37]\], using the above selected evolutionary substitution model and implementing two runs of four Markov Chain Monte Carlo (MCMC) processes with 7 million generations for *cytb* and 15 million generations for *GHI*, sampling every 100 generations in both cases.  The difference in chain length was because *GHI required* more generations to obtain sufficient values of effective sample size (\> 200).  We evaluated the convergence of the log-likelihood (-InL) values of the two runs, with 10% of reconstructions discarded as burn-in, to construct the consensus tree (σ \< 0.005).  We used the posterior probabilities obtained based on a confidence limit of 95% (highest posterior density-HPD) to evaluate the support of nodes.  We used a incongruence length difference test (Partition homogeneity test) \[[@CR38]\] in PAUP\* 4.0b10 \[[@CR39]\] to evaluate the significance of conflict among data sets.  We used 1000 resampling characters, and performed concatenated analyses of nuclear and mitochondrial data for ML and BI with the parameters previously mentioned.  We constructed haplotype networks for both genes using PopART (available at htt://popart.otago.ac.nz) and applying the Median-Joining method \[[@CR40]\]. 

Divergence times estimation and ancestral area reconstruction {#Sec6}
-------------------------------------------------------------

We estimated divergence times by a coalescent-based method using the software BEAST 1.8.0 \[[@CR41]\], based on the described evolutionary substitution model for both genes (see Additional file [1](#MOESM1){ref-type="media"}).  We used a relaxed molecular clock, using \*Beast (including nuclear an mitochondrial data) with Yule speciation tree model.  The clock was calibrated using the fossil record (in hard minimum bound for Lognormal distribution), which was located in the root of the *Pantosteus plebeius-nebuliferus* species-group.  This was represented by fossils of *Pantosteus asitus*, dated from the Miocene-Pliocene (7.5-2.5 Mya) \[[@CR42]\].  We conducted another molecular clock analysis using the mutation rate range for *cytb* estimated for teleosts of 0.76-2.2%/Myr with uniform distribution in BEAST \[[@CR43]--[@CR46]\]. 

We implemented three independent \*Beast analyses, each with 70 000 000 generations, sampling every 1000 generations.  After evaluating the posterior parameter values based on effective sample size, and according to the convergence data, using Tracer 1.6 \[[@CR47]\], we discarded 10% of runs and combined the three analyses to construct a maximum clade credibility tree using BEAST modules (LogCombiner and TreeAnnotator, respectively). 

We performed an ancestral area reconstruction with RASP 3.2, 20160719 \[[@CR48]\], applying Statistical-Dispersal Vicariance Analysis (S-DIVA) \[[@CR49], [@CR50]\] and the Statistical Dispersal-Extintion-Cladogenesis model (S-DEC) \[[@CR51]\].  In both cases we used the chronograms obtained in BEAST to resample every 3000 samples from the total set to obtain a subset of 70 000 trees.  We considered nine areas in the analysis corresponding mainly to the drainage basins sampled The areas included were the basins: A, Upper North Rio Grande (Palomas); B, Casas Grandes; C, Santa Maria; D, Del Carmen; E, Fuerte; F, Upper South Rio Grande (Conchos); G, Piaxtla; H, Nazas; and I, Mezquital.  We applied an unconstrained model allowing any combination of geographic range in adjacency matrix with a maximum of five adjacent hydrographic basins, configuring the dispersal rate with the same probability for each population (basin) of the *Pantosteus plebeius-nebuliferus* species-group. 

Species delimitation analysis {#Sec7}
-----------------------------

We performed species delimitation tests using Bayesian Phylogenetics and Phylogeography software BPP 3.4 \[[@CR52]\] to generate posterior probabilities for two species-delimitation models in the *Pantosteus plebeius-nebuliferus* species-group,.  We used the reversible-jump Markov chain Monte Carlo (rjMCMC) methodology \[[@CR53]\] to implement an species delimitation analysis *A10* in BPP using a fixed guide tree (in this case a species tree) for a species-delimitation model assuming four species according with the main clades observed, which was obtained in the divergence time estimation as mentioned above.  We provided an additional guide tree for a species-delimitation model assuming two species, where we assigned individuals to lineages defined by recognized species *P. plebeius* and *P. nebuliferus* according with the distribution described for these species \[[@CR16]\].  We performed an analysis *A11* for unguided species delimitation analysis to determine the best species-delimitation model compared to *a priori* models mentioned above under *A10* analyses. 

Including both molecular markers, we analyze three combinations of parameters for population size (*Θ*) and species divergence time (τ): In large population size with deep divergences (*Θ α*=1, *β*=10; τ *α*=1, *β*=10); small population size with low divergences (*Θ α*=2, *β*=2000; τ *α*=2, *β*=2000); and large population size with low divergences (*Θ α*=1, *β*=10; τ *α*=2, *β*=2000).  We ran the analyses for 500 000 generations, sampling each five generations, with a final burn-in of the first 25 000. 

Results {#Sec8}
=======

We obtained a total of 92 DNA sequences of the cytb gene (1140bp) (Genbank accession numbers: MG203619-MG203705) and 83 of the 3^rd^ intron of GHI (654bp) (Genbank accession numbers: MG203706-MG203869) for Pantosteus plebeius (*n* = 63 of cytb and *n* = 49 of GHI) and P. nebuliferus (*n* = 29 of cytb and *n* = 33 of GHI) from 22 localities in 8 drainage basins, representing most of the species distribution range \[[@CR16]\].  We sequenced 87and 82 sequences respectively for cytb and GHI and obtained5 of cytb and 1 for GHI from GenBank (Fig.  [1](#Fig1){ref-type="fig"} and Table [1](#Tab1){ref-type="table"}).  We did not detect any significant recombination in the nuclear GHI sequences (*p* \> 0.05). 

Genetic diversity {#Sec9}
-----------------

Although the haplotype diversity in all basins was relative high (h = 0.829-0.981) for *cytb* and *GHI*, the nucleotide diversity and the proportion of segregating sites were low for both markers (π = 0.001-0.010 and Θ~S~ = 0.0005-0.010) (Table [2](#Tab2){ref-type="table"}).  The populations with the lowest nucleotide diversity (π = 0.003) were those from Fuerte basin.  The Casas Grandes and Nazas basin populations showed the highest number of haplotypes.  However, the number of haplotypes was different in each marker by population as result of the difference in variation that each molecular marker has. Table 2Genetic diversity for *cytb*\|*GHI* markers in *Pantosteus plebeius-nebuliferus* populationsBasinsHnπhΘ~S~Nazas12\|180.010\|0.0030.850\|0.8290.009\|0.003Mezquital8\|30.008\|0.0010.924\|0.3420.010\|0.001Piaxtla1\|3\--\|0.008\--\|0.464\--\|0.001Fuerte7\|20.003\|0.0010.873\|0.4710.006\|0.0005Upper South R. Grande (Conchos)2\|20.003\|0.0011\|0.5000.003\|0.001Upper North R. Grande (Palomas)1\|1\-\-\-\-\--Del Carmen5\|20.005\|0.0010.800\|0.5030.006\|0.0005Santa Maria7\|70.005\|0.0060.964\|0.7920.005\|0.003Casas Grandes13\|60.006\|0.0020.981\|0.6830.008\|0.002*Hn* number of haplotype, π nucleotide diversity, *h* haplotype diverisity, Θ~*S*~ proportion of segregating sites

Genetic distances and structure {#Sec10}
-------------------------------

The highest absolute genetic distances (AGD) for *cytb* found between clades were between the Mezquital (Clade IV) and Piaxtla/Nazas basins (Clade III) (D~P~ = 9.2%) (Table [3](#Tab3){ref-type="table"}), and higher AGD than 5% between all clades was observed.  Based on *GHI* the AGD highest distances (D~P~ = 0.9%) were between the Mezquital and Fuerte and Del Carmen populations, and between the Piaxtla and Del Carmen and Fuerte populations.  Between all clades AGD distances were lower than one percent.  The lowest AGD were observed between populations of Casas Grandes and Santa Maria populations for *cytb* (D~P~ = 1.3%) and between the Palomas and Casas Grandes and Santa Maria for *GHI* (D~P~ = 0.1%), both comparisons within the Clade I (Guzman hydrographic system).Table 3Absolute pairwise un-corrected *p*-distances D~P~ for *cytb* (under the diagonal) and *GHI* (above the diagonal using pairwise-deletion) of populations of *Pantosteus plebeius-nebuliferus* species-groupMezquitalFuerteUSR Grande (Conchos)Casas GrandesSanta MariaDel CarmenPiaxtlaNazasUNR Grande (Palomas)Mezquital\--0.90.80.80.80.90.70.80.6Fuerte8.7\--0.10.70.70.80.90.70.5UNR Grande (Conchos)8.80.6\--0.60.60.70.80.70.5Casas Grandes8.65.75.9\--0.20.20.70.60.1Santa Maria8.16.66.62.6\--0.20.80.70.1Del Carmen7.95.75.71.32.0\--0.90.70.2Piaxtla8.97.67.27.67.97.5\--0.40.6Nazas9.38.17.77.87.97.81.4\--0.5USR Grande (Palomas)8.56.36.31.72.61.67.48.0*UNR* Upper North River, *USR* Upper South River

AMOVA revealed significant genetic structures among populations (Φ~ST~ = 0.9028 for *cytb*; Φ~ST~ = 0.6904 for *GHI*; *p* \< 0.05) (Table [4](#Tab4){ref-type="table"}).  The highest genetic structure with SAMOVA was found when K = 8 (Φ~CT~ = 0.8784 for *cytb*; Φ~CT~ = 0.7460 for *GHI*, *p* \< 0.05), and each basin was clustered in an independent group, with the exception of one population of the current Upper South Rio Grande, the Conchos population, which was clustered with the Fuerte population.  In the case of *GHI*, no significant genetic structure was observed when groups of populations (basins) recognized for *P. plebeius* and *P. nebuliferus* were compared (Φ~CT~ = 0.1707, *p* \> 0.05).Table 4Spatial Analysis of Molecular Variance (SAMOVA) of the *Pantosteus plebeius-nebuliferus* species-group baes on *cytb*\|*GHI*.  Only the results for K = 2, 4, and 8 are shown, as they are meaningful in phylogenetic, systematic, and biogeographic termsCorrespondenceKGroupsΦ~CT~Φ~SC~Φ~ST~One genetic pool (via AMOVA)\-\-\-\-\--0.9028\*\|0.6904\**Pantosteus plebeius*/*Pantosteus nebuliferus*2(Palomas/Casas Grandes/Santa Maria/Del Carmen/Conchos/Fuerte/Mezquital) (Piaxtla/Nazas)0.3545\*\|0.17070.8742\*\|0.6606\*0.9188\*\|0.7185\*Phylogenetic inference (Clades)4(I) (II) (III) (IV)0.7685\*\|0.6166\*0.6324\*\|0.2715\*0.9149\*\|0.7207\*Best genetic structure identified without geographic information8(Pamolas) (Casas Grandes) (Santa Maria) (Del Carmen) (Conchos/Fuerte) (Piaxtla) (Nazas) (Mezquital)0.8784\*\|0.7460\*0.2051\*\|0.21560.9033\*\|0.6912\**\*p* \< 0.05

High Φ~ST~ and significant pairwise values were observed among basins in both molecular markers, with Φ~ST~ ranging form 0.461 to 0.968 for *cytb* and form 0.050 to 0.907 for *GHI* (*p* \< 0.5) (Table [5](#Tab5){ref-type="table"}).  The highest value was observed between the populations of Piaxtla and Fuerte for *cytb*, while between Piaxtla and Del Carmen populations for *GHI.*  The lowest but significant Φ~ST~ was between Piaxtla and Nazas for *cytb* and between Del Carmen and Santa Maria for *GHI*.  In contrast to *cytb* Φ~ST~ values, no significant Φ~ST~ between populations of Santa Maria and Casas Grandes were observed for *GHI*.Table 5Genetic differentiation using pairwise Φ~ST~ for *cytb* (under the diagonal) and *GHI* (above the diagonal) among basinsNazasPiaxtlaMezquitalFuerteConchosCasas GrandesDel CarmenSanta MariaPalomasNazas\--**0.5320.7420.722**0.656**0.6400.7230.572**0.619Piaxtla**0.461**\--**0.8750.897**0.914**0.7830.9070.596**0.900Mezquital**0.8960.931**\--**0.896**0.899**0.8010.9030.657**0.890Fuerte**0.8990.9680.936**\--0.050**0.7770.8870.635**0.868Conchos0.8740.9900.9170.473\--0.7100.8940.4121.000Casas Grandes**0.8890.9400.9230.921**0.911\--**0.314**0.0130.327Del Carmen**0.8880.9600.9220.938**0.934**0.649**\--**0.178**0.789Santa Maria**0.8840.9510.9140.936**0.921**0.7950.777**\--0.200Palomas0.8721.0000.9080.9510.9560.6550.7630.780\--Bold format = *p* \< 0.05 after Bonferroni correction

Phylogenetic relationships and haplotype network analysis {#Sec11}
---------------------------------------------------------

The incongruence length difference test did not show significant differences (*p* \< 0.05), indicating that both molecular markers presented the same phylogenetic signal.  Bayesian and Maximum Likelihood phylogenetic analyses for the *cytb* and *GHI* sequences (see Additional files [2](#MOESM2){ref-type="media"} and [3](#MOESM3){ref-type="media"}) and concatenated analysis (Fig.  [2](#Fig2){ref-type="fig"}) recovered the same topology, with four highly supported clades (main clades).  Clade I included samples from the Guzman hydrographic system (Casas Grandes, Santa Maria, and Del Carmen basins) and the Upper North Rio Grande system (Palomas).  Clade II comprised the populations of the Fuerte basin and Conchos River (Upper South Rio Grande basin).  Clade III, corresponding to the Nazas and Piaxtla basin populations and Clade IV consisted of the population from the Mezquital basin. Fig.  2Phylogenetic reconstruction *Pantosteus plebeius-nebuliferus* species-group based on the concatenated data of *cytb* gene and *GHI* region.  The numerical values in nodes represent the posterior probabilities and the bootstrap values for Bayesian Inference and Maximum Likelihood, respectively

The same reciprocally monophyletic haplogroups were recovered in the haplotype networks reconstruction for both molecular markers, corresponding to the four main clades obtained by phylogenetic analysis, but differ in the relationships.  In *cytb*, the Mezquital population (Clade IV) is related (72 mutations steps or MS) to the Guzman System/Palomas populations (Clade I), whereas the Nazas/Piaxtla populations (Clade III) are more closely related to the Conchos/Fuerte populations (72 MS) (Clade II) than to the Mezquital population (Fig.  [3a](#Fig3){ref-type="fig"} and Fig.  [3b](#Fig3){ref-type="fig"}).  The number of MS separating the population of Fuerte from the Conchos River; and the population of Nazas from Piaxtla River, was 3 MS in both cases.  We did not include as MS the median vectors or hypothetical haplotypes implied by the haplotype network reconstruction. Fig.  3Haplotype network for *Pantosteus plebeius-nebuliferus* species-group based on *cytb* (**a**) and *GHI* (**b**).  The numbers inside the circles represent the number of mutations steps between haplotypes.  The grey circles represent median vector or hypothesised sequences which is required to connect existing sequences within the network

In case of *GHI*, the range of MS between clades was 1-2 MS, geographic correspondence was observed, and shared haplotypes were identified between populations that are closely related phylogenetically (Casas Grandes/Del Carmen/Santa Maria/Palomas, Fuerte/Conchos and Piaxtla/Nazas) (Fig.  [3b](#Fig3){ref-type="fig"}). 

Divergence times and ancestral area reconstruction {#Sec12}
--------------------------------------------------

The global likelihood scores with the highest probabilities were obtained using the S-DEC approach (Fig.  [4](#Fig4){ref-type="fig"}).  However total congruence in ancestral areas predicted was observed between both analyses.  The scenario showed by the two models (S-DIVA and S-DEC) inferred an ancestral area (ABCDF) of the most recent common ancestor (MRCA) of the *Pantosteus plebeius-nebuliferus* species-group formed by the Guzman system, and the Upper South and North Rio Grande (Conchos and Palomas, respectively).  From this first vicariant event dated at *ca.*  5.1-3.9 Mya (HPD ≥ 95%) to the upper Pliocene, a succession of 7 vicariant and 4 dispersal events were identified during the historical distribution of the species (Fig.  [4](#Fig4){ref-type="fig"}).Fig.  4Divergence times estimation (based on species tree), S-DIVA and S-DEC analyses of the *Pantosteus plebeius-nebuliferus* species-group.  Red labels represent the divergence times in millions years ago.  The bold names represent the ancestral area estimated.  Dispersal events and vicariant events are represented with - and \| symbols, respectively.  The \* represent the value for posterior probability of 100% for each clade tested in the most sensitive species-delimitation model (four species) with values of *Θ α*=1, *β*=10; *τ α*=1, *β*=10 in the species delimitation analysis *A10* and *A11*

Species limits {#Sec13}
--------------

The species-delimitation model assuming four species was strongly sensitive to low parameter of *Θ* and τ, showing high posterior probability for the nodes in *A10* and was the species-delimitation model in *A11* analysis with the highest probability (P=1 S=4) (Fig.  [4](#Fig4){ref-type="fig"}).  The four species supported correspond to the main clades identified in this work for *Pantosteus plebeius-nebuliferus* (Fig.  [2](#Fig2){ref-type="fig"}).  The model of species-delimitation reducing the number of species to two (*P. plebeius* and *P. nebuliferus*) showed low posterior probabilities for the three combinations of parameters of *Θ* and τ (results not showed). 

Discussion {#Sec14}
==========

The main forces shaping cladogenetic events in the *Pantosteus plebeius-nebuliferus* species-group are consistent with geological processes occurring since the Pliocene associated with the formation and evolution of the SMOC, central North Mexico and the tectonic activity of the Rio Grande Rift.  Our results revealed that the species-group comprises more than two independent evolutionary lineages with significant genetic structures between them in nuclear and mitochondrial markers, indicating a long history of isolation.  Within clades mitochondrial information showed high divergence and several mutational steps between clades, unlike the nuclear marker that showed shared haplotypes between some basins; and low number of mutational steps between clades, reflecting the combined effects of low variation and incomplete lineage sorting on the nuclear loci \[[@CR54]\].  As has occurred with the genus *Moxotoma* \[[@CR12]\] and *Campostoma ornatum* \[[@CR3]\], relevant vicariant events and probably river captures may have played an important role in range expansion shaping the spatial distribution of genetic variation in some drainage basins.  Genetic variation within drainage basins seems to have been influenced by ecological characteristics of the species and climatic fluctuations since the Pleistocene.  Further studies should seek to identify the roles of historical, intrinsic, and anthropogenic influences on genetic differentiation within basins \[[@CR55]\]. 

Historical biogeography of pantosteus plebeius-nebuliferus {#Sec15}
----------------------------------------------------------

### Pliocene events {#Sec16}

We dated the most recent common ancestor of all clades of the *Pantosteus plebeius-nebuliferus* species-group to the Pliocene (5.1-3.9 Mya, HPD ≥ 95%) (Fig.  [5](#Fig5){ref-type="fig"}), in an ancestral area comprised by the proto-Guzman system, proto-Upper South Rio Grande (Conchos) and proto-Upper North Rio Grande (Palomas).  We identified a range expansion of the species-group to the Fuerte basin, which probably resulted from a basin catchment, allowing faunal interchange as proposed for *Gila pulchra* \[[@CR10]\], and *Codoma ornata* \[[@CR11]\].  These events occurred during the tecto-volcanic episodes in SMOC evolution, including repeated volcanic activity called "alkaline basalt events" \[[@CR4], [@CR56], [@CR57]\].  Another dispersal event from the ancestral areas to Nazas was proposed.  These events could have been promoted by alkaline basalts events during the Mapimi Bolson formation (Fig.  [4](#Fig4){ref-type="fig"} and Fig.  [5](#Fig5){ref-type="fig"}) \[[@CR12], [@CR58], [@CR59]\], which is partially congruent with patterns established for *Codoma ornata* \[[@CR11]\] and *Campostoma ornatum* \[[@CR7]\].Fig.  5Biogeographical scenario for the *Pantosteus plebeius-nebuliferus* species-group in the North Mexico

These dispersal events were followed by a vicariant event that isolated the ancestor of the Proto Guzman system/Palomas (Clade I) *ca.*  4.5 Mya, associated with the volcanism in the Babicora-Bustillos sector, disrupting the connection between the proto-Fuerte/Conchos and the proto-Guzman systems/Palomas by the periodic accumulation of sediments in the Conchos Valley \[[@CR12]\].  This ancient isolation event is reflected in the high genetic divergence between clades I and II in *cytb* (D~P~ = 6.3-7.0%, Φ~ST~ = 0.911-0.938 *p* \< 0.05; Table [3](#Tab3){ref-type="table"}) and in the well-structured haplogroups for both markers (Fig.  [3](#Fig3){ref-type="fig"}).  The split between clade I with respect to the clades II, III, and IV of the *P. plebeius-nebuliferus* species-group, is geographically, genetically and in some cases temporarily consistent with other co-distributed species complexes of fishes such as *Cyprinella* spp.  \[[@CR60]\], *Campostoma ornatum* \[[@CR3], [@CR61]\], *Cyprinodon* spp*.*  \[[@CR6], [@CR62]\], *Gila* spp.  \[[@CR10]\], *Codoma ornata* \[[@CR11]\] and *Moxostoma* cf. * austrinum* from Conchos \[[@CR12]\].  This finding also corroborates the divergence time of the Rio Grande *Pantosteus* spp.  population with respect to the Nazas population in a previous work \[[@CR63]\]. 

The Mapimi Bolson desertification and the alkaline basalt activity at the central part of the SMOC, was likely responsible for the disconnection that promoted the vicariant event of the ancestor of the proto-Nazas/Piaxtla (Clade III)/Mezquital basin populations (Clade IV) from the proto-Fuerte/Conchos populations in the Pliocene, an event that is consistent with results found for *Codoma ornata* results \[[@CR11]\].  Two range expansion events of the species-group from the Nazas to Piaxtla and Mezquital basins following the vicariant event, were found in the same period, which could be the result of river catchments associated with alkaline basalt events in the SMOC \[[@CR4], [@CR12], [@CR57]\], an event that is consistent with results found in other aquatic fauna, such as *Campostoma ornatum* \[[@CR3]\] and semi-aquatic snake species of the genus *Thamnophis* \[[@CR64]\]. 

Due to the physiographic complexity and the occurrence of several overlapping geologic events in the region of Nazas, Piaxtla and Mezquital, is difficult to identify specific geologic events involved in the cladogenesis of clades III and IV occurred *ca.*  2.6 Mya.  The isolation of the proto-Nazas/Piaxtla from the Mezquital populations possibly had its source in the alkaline volcanic activity occurring in the region during the Pliocene, especially south of Laguna Santiaguillo in the Nazas basin, in the southern region of the SMOC \[[@CR4]\] (Fig.  [4](#Fig4){ref-type="fig"}) separating the tributary El Tunal (Mezquital Basin) from the Nazas basin occurring in the Late Miocene and Pliocene \[[@CR2], [@CR65], [@CR66]\]. 

### Pleistocene isolation {#Sec17}

The origin and evolution of the three basins within the Guzman system date from the Late Cretaceous and Cenozoic to the Late Pleistocene, and its complexity could explain the species diversity and genetic divergence of fishes in the region.  During the Early Pleistocene, the activity of the Rio Grande Rift and the glacial stage (Kansan, Middle Pleistocene) \[[@CR67]\] were likely involved in the fragmentation of the Cabeza de Vaca paleolake and the formation of the Palomas pluvial paleolake, a relict of the Cabeza de Vaca that was fragmented into several endorheic rivers, lagoons, and springs.  Among these, the Casas Grandes and Upper Rio Grande basins were firstly isolated from Del Carmen and Santa Maria (*ca.*  1.8 Mya) by the formation of the Los Muertos Bolson \[[@CR60], [@CR68]--[@CR70]\] (Fig.  [5](#Fig5){ref-type="fig"}).  The isolation of the Upper North Rio Grande, including its flow to the Gulf of Mexico and capture of the Conchos River, was dated to the Early Pleistocene \[[@CR71]--[@CR74]\], but the tectonic subsidence in the Rio Grande Rift and desertification of Los Muertos Bolson has continued during the Middle Pleistocene and the present \[[@CR75]--[@CR77]\] allowing the fragmentation of the Guzman system.  This paleohydrological pattern is reflected by our data, since we date the isolation of Palomas from Casas Grandes as *ca.*  1.1 Mya and 1.2 Mya for the isolation of Del Carmen from Santa Maria (Fig.  [4](#Fig4){ref-type="fig"}).  This is also consistent with previous studies that found low genetic divergences between populations of *Pantosteus plebeius* from the Mimbres River (type locality of the species) and *P. plebeius* from the Palomas River (Upper North Rio Grande) \[[@CR22]\]. 

The three basins within the Guzman system were isolated and reconnected several times during the Pre-Illinoian, Illinoian, and Pre-Holocene periods \[[@CR78]--[@CR80]\].  Remixing gene pools and low genetic structures also were observed in other co-ocurring fishes along the Guzman basin, including *Gila* spp.  \[[@CR10], [@CR81]\] and *Cyprinodon* spp.  \[[@CR6], [@CR62]\].  However, the high interbasin genetic distances and the lack of shared haplotypes for *cytb* is evidence of the ancestral isolation that occurred in Del Carmen and Casas Grandes populations of *Campostoma ornatum* (D~P~ = 2.9%; Φ~ST~ = 0.935) \[[@CR3]\], and the significant genetic differentiation (F~ST~ = 0.380) in *Cyprinella formosa* \[[@CR60], [@CR69], [@CR82]--[@CR84]\] which represents an isolation event that also occurred in the Pleistocene.  In contrast, the low genetic differentiation and haplotypes shared among basins for *GHI* are evidence for the low variation and of the incomplete lineage sorting or a later genetic mixing between these basins.  More studies are necessary to explain the configuration of the ichthyofauna of the Guzman system. 

Finally, the vicariant events involved in the isolation of the Nazas from the Piaxtla and the Conchos from the Fuerte populations, were estimated to have occurred 1.4 Mya and 0.6 Mya in the Pleistocene (Fig.  [4](#Fig4){ref-type="fig"}), respectively, and are associated with the most recent volcanic eruptions and tectonic movements in the SMOC \[[@CR85]\].  This supports the results obtained for *Campostoma ornatum* \[[@CR3]\], *Gila pulchra* (Conchos River & Fuerte basin) \[[@CR10]\], and *Codoma ornata* (upper Conchos River & Fuerte basin) \[[@CR11]\] that show low genetic differentiation as well as shared haplotypes among basins as we observed in the nuclear marker (Fig.  [3](#Fig3){ref-type="fig"}). 

### Within-basin differentiation {#Sec18}

A high number of mutational steps and genetic intrabasin distances where estimated for *cytb* in the *Pantosteus plebeius-nebuliferus* species-group*.*  Surprisingly, moderate between-haplotype MS values (MS = 4-7) were found within the Mezquital, Nazas, Fuerte, Casas Grandes, and Santa Maria basin populations than between the Nazas and Piaxtla (MS = 3) and Conchos and Fuerte populations (MS = 3), reflecting considerable genetic differentiation in *cytb* (Fig.  [3](#Fig3){ref-type="fig"} and Tables [3](#Tab3){ref-type="table"}, [4](#Tab4){ref-type="table"}, [5](#Tab5){ref-type="table"}).  The intrapopulation genetic structure in those drainage basins could be explained in part by life history traits.  Migration episodes associated with environmental fluctuations on spawning sites and food availability have been extensively described in other species of Catostomidae \[[@CR86]--[@CR93]\].  Accordingly, the historic climate fluctuations in the area \[[@CR94]--[@CR96]\] could promote similar fragmentation episodes in *P. plebeius-nebuliferus* populations, which could be responsible for the isolation of migrants in separate regions within the basins. 

Taxonomic implications and conservation considerations {#Sec19}
------------------------------------------------------

The *Pantosteus plebeius-nebuliferus* species-group comprises a set of distinct genetic groups, suggesting that it is composed of at least four clades with unique and independent evolutionary histories.  We found cladogenesis to be associated with the high tecto-volcanic activity and climate fluctuations from the Pliocene to the present in North Mexico.  This is supported by the phylogenetic relationships, species tree, species-delimitation test including both molecular markers.  Also is supported by the genetic divergences observed between clades for *cytb* (see Additional file [2](#MOESM2){ref-type="media"}), which are higher than those reported among sibling vertebrate groups for *cytb* (\~2%) \[[@CR97]--[@CR100]\] and between sister species within the family \[[@CR98]\].  In the case of *GHI*, the genetic distances estimated between species of Catostomidae are 3.3% \[[@CR101]\] taking into account different genera of the family, and which were higher than estimated in this work (Table [3](#Tab3){ref-type="table"}, D~P~\<1%).  However, although the nuclear *GHI* showed a low number of mutational steps and genetic distances (Fig.  [3](#Fig3){ref-type="fig"} and Table [3](#Tab3){ref-type="table"}), we detected significant structure (Table [4](#Tab4){ref-type="table"}) and recovered the four clades in the phylogenetic inference (see Additional file [3](#MOESM3){ref-type="media"}) supporting the mitochondrial results.  The genetic structure and divergences in the species-group were similar to those observed in other fish taxa that partially or completely occur in basins of northwest Mexico and the southwestern USA, such as *Cyprinodon* spp.  \[[@CR6], [@CR62]\], *Gila s*pp.  \[[@CR10], [@CR81]\], *Cyprinella* spp.  \[[@CR60], [@CR69], [@CR84]\], or even among populations of other considered species complexes as *Campostoma* spp.  \[[@CR3], [@CR61]\], and *Moxotoma* spp.  \[[@CR12]\].  Thus, we recommend an integrative taxonomic revision of the *Pantosteus plebeius-nebuliferus* to recognize the four independent evolutionary lineages as different species. 

The population of the Piaxtla basin was recognized as *Pantosteus plebeius* \[[@CR16]\] and suggested as an independent evolutionary lineage \[[@CR20]\], our results show that this population is closely related to *P. nebuliferus*, as was also demonstrated by \[[@CR63]\].  Low genetic divergences between populations of *P. plebeius* from the type locality Mimbres River and from the Upper North Rio Grande have been observed by other authors \[[@CR22], [@CR102]\].  Our results show that Paloma River is closely related with Guzman system populations, thus we suggest that Clade I must be considered populations of the recognized species *P. plebeius*.  The Clade II (Fuerte basin and Conchos River populations) and the Clade IV (Mezquital River population) do not represent populations of any species recognized of the species-group studied, thus we suggest that the two clades could represent two independent evolutionary lineages that must be described and recognized as new species of the genus *Pantosteus*. 

The identification of the Evolutionary Significant Units (ESU) is of particular interest especially in species considered threatened \[[@CR103]\].  Such is the case of *Pantosteus plebeius-nebuliferus* species group, which is listed as in danger of extinction in the Norm-059-SEMARNAT-2010 \[[@CR104]\].  We identified the group of populations whitin independent evolutionary lineages, which warrant separate management or priority for conservation.  Based on the high genetic divergences; number of mutation steps; phylogenetic relationships; shared haplotypes; reproductive isolation and the subsequent absent exchangeability of populations, we suggest nine ESUs corresponding to each basin where the *P. plebeius-nebuliferus* species-group inhabit: (1) Upper North Grande River, (2) Casas Grandes, (3) Santa Maria, (4) Del Carmen, (5) Fuerte, (6) Upper South Grande River basin (Conchos River), (7) Piaxtla, (8) Nazas and (9) Mezquital.  The definition of these nine ESUs should aid in establishing conservation measures of the species-group and of the evolutionary lineages of which they are configured. 

Conclusions {#Sec20}
===========

Our study supports at least four independent evolutionary lineages in the *Pantosteus plebeius-nebuliferus* species-group.  Their biogeographic history is linked to geological events occurring since the Pliocene, associated with the formation and evolution of the SMOC, central North Mexico and the tectonic activity of the Rio Grande Rift.  Other freshwater fishes already studied and occurring partially or completely in the same range of northwest Mexico and the southwestern USA have similar genetic structures and divergences.  We recommend an integrative taxonomic revision of the *Pantosteus plebeius-nebuliferus* species-group and propose nine ESUs for conservation purposes which are restricted to individual basins where the species-group is distributed.  As an endangered species, we hope our survey contributes to the management and preservation of these fishes. 

Additional files
================

 {#Sec21}

Additional file 1:Evolutionary substitution model and estimated parameters for *cytb* and *GHI* by Akaike Information Criterion. ( DOC 35 kb) Additional file 2:Phylogenetic tree (Bayesian Inference/Maximum Likelihood) based on mitochondrial gene *cytb* for the *Pantosteus plebeius-nebuliferus* species-group. ( EPS 1046 kb) Additional file 3:Phylogenetic tree (Bayesian Inference/Maximum Likelihood) based on nuclear gene *GHI* for the *Pantosteus plebeius-nebuliferus* species-group. ( EPS 1151 kb)

AGD

:   Absolute Genetic Distances

AMOVA

:   Analysis of Molecular Variance

BI

:   Bayesian Inference

*ca*

:   Approximately

CPUM

:   Fish Collection of the Universidad Michoacana de San Nicolás de Hidalgo

*cytb*

:   Cytochrome b

D~P~

:   Genetic *p*-distances

ESU

:   Evolutionary Significance Unit

GHI

:   3rd Intron of Growth Hormone copy I

GTR+G

:   General Time-Reversible + Gamma

h

:   Haplotype diversity

HN

:   Number of haplotype

HPD

:   Highest Posterior Density

-lnL

:   Log-likelihood

MCMC

:   Markov Chain Monte Carlo

ML

:   Maximum Likelihood

MRCA

:   Most Common Ancestor

MS

:   Mutation Steps

Mya

:   Million years ago

Myr

:   Million years

PCR

:   Polymerase Chain Reaction

SAMOVA

:   Spatial Analysis of Molecular Variance

S-DEC

:   Statistical Dispersal-Extintion-Cladogenesis

S-DIVA

:   Statistical-Dispersal Vicariance Analysis

SMOC

:   Sierra Madre Occidental

TrN+G

:   Tamura-Nei + Gamma

UNR

:   Upper North River

USA

:   United States of America

USR

:   Upper South River

*Θ*

:   Population size

Θ~S~

:   Proportion of segregating sites

π

:   Nucleotide diversity

τ

:   Species divergence time

Φ~ST~

:   Pairwise fixation indices

We thank Rogelio Rosas Valdez and Hugo Mejia Madrid for field support in obtaining fish samples.  Thanks to Javier Pérez Tris and Silvia Perea for valuable comments, and to the reviewer Tom Jameson from Chester Zoo and reviewers from The Lucidus Consultancy, Alan Pike and Kathleen Hills, for English corrections and suggestions. 

Funding {#FPar1}
=======

This work was partial supported by the Coordinación de la Investigación Científica (CIC-UMSNH 2014), the Ministerio de Economía y Competitividad of Spain (grant reference CGL2013-41375-P).  Both resources were used to finance fieldwork, laboratory (DNA isolation and amplification of molecular markers) and sequencing service, as well as for the publication of this work.  The Consejo Nacional de Ciencia y Tecnología CONACYT with sabbatical grant and the scholarship to ODD (grant number: 260596) and DKCS (number grant: 254270/411621) respectively. 

Availability of data and materials {#FPar2}
==================================

The datasets generated and/or analyzed during the current study are included in the research article (and its additional files).  Also, all DNA alignments and phylogenetic trees are available via Treebase (<http://purl.org/phylo/treebase/phylows/study/TB2:S23572?x-access-code=31a6674c7214877743d545e8ded23145&format=html>).  All DNA sequences obtained in this work are available on GenBank under the Accession Numbers MG203619-MG203705 for *cytb* gene, and MG203706-MG203869 for *GHI* region. 

All authors contributed intellectually to the work and approved the submission.  ODD, RRP, LTM, YHD, ID and DKCS participated in fieldwork.  DKCS, LTM and JRPB carried out the work lab.  DKCS, LTM, JRPB, ODD and ID analyzed the data and wrote the manuscript.  All authors have read and approved the manuscript for submission. 

Ethics approval and consent to participate {#FPar3}
==========================================

This research has been performed within an appropiate ethical frameworks.  Field and laboratory protocols used in this study, including sampling procedures, were reviewed and approved by the Mexican Ministry of Environmental and Natural Resources (SEMARNAT), under collection permit number PPF/DGOPA-262/17.  Approval by an ethics committee was not necessary because this research did not include animal experimentation, and the samples were obtained using tricaine mezylate (MS-222) in fishes as anesthetic, according with the Official Mexican Norm NOM-051-ZOO-1995 and NOM-033-SAG/ZOO-2014 for humanitarian treatment in the mobilization of animals. 

Consent for publication {#FPar4}
=======================

Not applicable. 

Competing interests {#FPar5}
===================

The authors declare that they have no competing interests. 

Publisher's Note {#FPar6}
================

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.