Samples were collected through field work and from specimens in captivity, with declared geographical origin. Furthermore, we included sequences from the GenBank database with the aim to incorporate northernmost populations (Southwestern USA) as well as C. totonacus. The list of sequences used in this work can be found in Suppl. material 1: Table S1 (GenBank accession numbers MN327706–MN327962, sequences obtained in this work) and S3(Crotalus molossus from Anderson and Geenbaum 2012 and outgroups). We sampled 84 specimens collected throughout their range from the Sierra Madre del Sur in Oaxaca, across central and northern Mexico, to the southern United States of America (Fig. 2a): 9 samples of C. m. molossus, 41 of C. m. nigrescens, 5 of C. m. oaxacus, 20 of C. basiliscus, and 8 of C. ornatus. The individuals were ethically handled, and tissue was obtained by ventral-vein blood extraction or ventral scales clipping (collecting permit: SGPA/DGVS/011587/17 SEMARNAT). All individuals were released at the point of capture after sampling. We provide collection locations of the tissue samples in Suppl. material 1: Table S1.

We extracted DNA from ventral scales and blood using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, California), following the manufacturer’s instructions. Ventral scales were additionally incubated in 60 µL DTT (DL-Dithiothreitol 5%) and 3X of proteinase K during lysis step. We amplified and sequenced mitochondrial genes (mtDNA), ATPase subunit 6 (ATPase 6), NADH dehydrogenase subunit 4 (ND4), cytochrome subunit b (cytb); as well as the nuclear recombination gene that activates protein 1 (RAG-1) (Suppl. material 1: Table S2). Those genes were chosen due to their levels of molecular polymorphism, which had been shown previously to be informative in inferring phylogenetic relationships in viperids and especially within rattlesnakes (Bryson et al. 2011; Anderson and Greenbaum 2012; Reyes-Velasco et al. 2013; Blair and Sánchez-Ramírez 2016; Myers et al. 2017; Meik et al. 2018).

We visualized the relative DNA quantity and quality on 1% agarose gels stained with the fluorescent stain GelGreen (Bradbury et al. 2010) under ultraviolet light. Amplification of the genes by Polymerase Chain Reaction (PCR) was carried out in a total volume of 13.05 µl containing 9.25 µl of ddH₂O, 1.25 µl of buffer (500 mM of KCl, 100 mM of Tris-HCl pH 9.1 at 20 °C and 0.1% TritonX-100), 0.4 µl of dNTP’s (2 mM), 0.5µl of forward primer and reverse primer (10 µM), 0.15 µl of Taq DNA polymerase (Vivantis) (5 U/µl) and 1 µl of 25–50 ng of DNA, with the following program in the Axygen Thermocycler: an initial denaturation of 94 °C for 5 minutes, followed by 35 cycles of 94 °C for 1 minute, 50 °C for 1 minute, 72 °C for 1 minute and a final extension of 72 °C for 10 minutes. Subsequently, 2µl of each sample was corroborated by electrophoresis comparing the bands with Axygen 100pb Ladder DNA (Corninig, NY). Purification and sequencing of samples was carried out at Laboratory of Genomic Services (LANGEBIO–CINVESTAV), Irapuato, Guanajuato, Mexico. Multiple independent PCR amplifications for random samples were also conducted and samples were sequenced twice to ensure reproducibility and correct readings.

We edited and cleaned sequences using CHROMAS 2.6.2 (Technelysium Pty Ltd.) and performed multiple sequence alignment with the MUSCLE algorithm (Edgar 2004) in MEGA 7.0.21 software (Kumar et al. 2016). Sequences were compared individually against the GenBank database, using the Basic Local Alignment Search Tool (BLAST) at the NCBI website (http://blast.ncbi.nlm.nih.gov/Blast.cgi), to ensure correct identification of amplifications. Subsequently, the sequences were traduced to proteins to inspect the presence of unexpected stop codons.

With the aim to consider the most complete taxonomic representation within black-tailed rattlesnakes, we developed two main molecular data sets: A) one consisting only of cytb, ND4 and ATPase6 loci (3-mtDNA dataset) and containing all six taxa (C. basiliscus, C. m. molossus, C. m. nigrescens, C. m. oaxacus, C. ornatus, and C. totonacus) and b) including the concatenated four loci, but containing missing data from the RAG1 (only available for 40 individuals, 15% missing data-dataset).

After considering different partitions’ schemes and datasets (Suppl. material 1: Fig. S3), we used the Akaike information criterion (AIC) from JMODELTEST 2.1.9 (Darriba et al. 2012) to select the best fitting model of sequence evolution for both concatenated datasets for the following phylogenetic analyses. Maximum likelihood was performed with PHYML 3.1 (Guindon et al. 2010), which is based on NNI+SPR for branch lengths and topology optimization. Clade support was assessed with 1,000 non-parametric bootstrap replicates. Bayesian inference was done with MRBAYES 3.2 (Ronquist and Huelsenbeck 2012) using default settings, and 4 chains were sampled every 1,000 generations for 50 million generations; convergence and stationarity within chains were visualized with TRACER 1.7 (Rambaut et al. 2018) and 25% of the generations were discarded as burn-in. We estimated the 50% majority-rule consensus topology and posterior probabilities for each node. In addition, we constructed a rooted haplotype network using the Neighbor-Net algorithm with SPLITSTREE 4.14.2 (Huson and Bryant 2006) based on the patristic distance corrected by the GTR+I+G model of evolution. Sequences from the genera Sistrurus (S. miliarius (Linnaeus, 1766) and S. catenatus (Rafinesque, 1818)), Agkistrodon (A. contortrix (Linnaeus, 1766) and A. piscivorus (Lacépéde, 1789)), and Crotalus (C. atrox (Baird & Girard, 1853), C. ruber (Cope, 1892), C. durissus (Linnaeus, 1758), C. tzabcan (Klauber, 1952) and C. simus (Latreille, 1801)) were used as outgroups in all phylogenetic analyses (Parkinson et al. 2000; Douglas 2006; Guiher and Burbrink 2008; Suppl. material 1: Table S3).

To investigate genealogical relationships considering solely the nuclear RAG1 locus, we inferred the phases of heterozygous genotypes with PHASE 2.1.1 (Stephens and Donnelly 2003), which reconstruct the most probable pair of alleles for each heterozygous individual, based on the frequencies of ambiguous haplotypes. Those phased alleles were used to construct an unrooted network with the Minimum-Spanning algorithm implemented with POPART 1.7 (Leigh and Bryant 2015). Since no C. totonacus individual was amplified for this marker, only five main lineages are represented following this scheme.

We performed Tajima’s D (Tajima 1989), Fu and Li’s F (Fu and Li 1993), and Fu’s FS (Fu 1997) tests to evaluate if the sequences within each dataset conformed to a neutral model of evolution in DnaSP 5.10 (Librado and Rozas 2009). We calculated three genetic diversity estimates for each dataset: the nucleotide diversity (p), haplotype diversity (h), and average number of differences among sequences (k). In order to estimate the genetic divergence among main clades, we estimated the average number of nucleotide substitutions per site between pairs of clades (Dxy, Nei 1987) and the number of net nucleotide substitutions per site between pairs of clades (Da, Nei 1987), with DNASP 5.10, considering the 15% missing data-dataset.

We used BEAST 1.8.4 (Drummond et al. 2012) to estimate the divergence times for the main lineages inferred. We applied a GTR+I+G model of evolution across all gene and codon positions, a Yule process tree prior, an uncorrelated lognormal relaxed molecular clock, and 50 million generations sampled every 1000^th generation, with 10% of the initial samples discarded as burn-in. Convergence, stationarity and an adequate effective sample sizes (>200 for each parameter) were visualized with TRACER 1.6 (Drummond et al. 2014). We considered two fossil-based references for set calibration points to carry out the phylogeny; for this, we included the most ancient fossil records for genus Sistrurus (late Miocene-Claredonian ~9 mya; Parmley and Holman 2007), specified with a normal logarithmic prior distribution of 9 (values in real space), mean of 0.01 and standard deviation of 1.0. (95% CI: 9.2–14.2 Mya). The second point is related to the most ancient fossil record for the genus Agkistrodon (Miocene-Late Hemphillian ~10 mya; Holman 2000); we specified a normal logarithmic previous distribution with displacement of 6 (values in real space), mean of 0.01 and standard deviation of 0.95 (95% CI: 6.2–10.8 Mya). Finally, we used calibration points to date the divergence between C. atrox and C. ruber and we specified a normal logarithmic distribution with displacement of 2.5 (normal in real space), mean of 0.01 and standard deviation of 0.50 (95% CI: 2.94–4.79 Mya). We ran several trial runs using the calibration points separately and combined, to ensure consistency among our results.

A Bayesian Skyride reconstruction was performed using BEAST 1.8.4. Due to the restricted number of unique haplotypes within the main clades (see the results below), we merged sequences from sister clades, as follows: C. basiliscus + C. m. molossus and C. ornatus + C. totonacus. We kept C. m. nigrescens apart due to its higher number of samples with respect the rest of the lineages. For C. m. oaxacus, despite its restricted number of unique haplotypes (five), we performed this analysis separately due to its basal position and divergence with respect to the rest of the clades. We sampled every 1000 generations along 10⁷ generations with 10% burn-in and using coalescent intervals (m) of 10. TRACER 1.7 was used to produce the skyride plot. Particularly, for C. m. oaxacus, we sampled every 1000 generations along 10¹⁰ generations, discarding 10% as burn-in and with coalescent intervals (m) set to 5.

To infer the ancestral distribution as well as the putative geographical origin of the black-tailed rattlesnakes we used RASP 4 (Yu et al. 2020) software, we performed the model selection using the BIOGEOBEARS method (Matzke 2013), which considers likelihood values and AIC to choose the best-fit model. This method implements DEC (Ree and Smith 2008), DIVA (Ronquist 1997), and BAYAREA (Landis et al. 2013) models, allowing for consideration of the jump dispersal (J parameter) (Matzke 2013). A total of six models (with and without J parameter) were compared using the summarized maximum clade credibility tree from BEAST analyses. To define the distribution areas required for the reconstruction of ancestral areas, 15 North American Terrestrial Ecoregions-Level III were considered (Wilken et al. 2011; Suppl. material 1: Table S4).

The 3-mtDNA based dataset successfully covered a total of 1,908 bp (ATPase 6=605, ND4=618, and cytb=685) from 60 samples of black-tailed rattlesnakes. This dataset yielded 493 polymorphic sites, 422 parsimony informative sites, and 58 haplotypes, yielding high haplotype diversity (h=0.999), moderate nucleotide diversity (π=0.074) and an average number of nucleotide differences between haplotypes (k) of 94.075. Tajima’s D (D= –0.3704; P>0.1), Fu and Li’s F (F= –0.17; P>0.1) and Fu’s F_S tests (F_S= –9.050; P>0.1) indicated no departures from neutrality. The model of sequence evolution selected for this dataset partition was GTR+I+G (pinv=0.339, α=0.86) and nucleotide frequencies were 32.91% A, 32.14% C, 9.83% G and 25.12% T. The nuclear gene (RAG-1=912 bp) yielded 182 polymorphic sites, 71 of them resulted as parsimoniously informative. This subsample (51 individuals) carried 30 phased haplotypes (h=0.893) with low nucleotide diversity (π=0.011). This molecular marker did not have enough molecular resolution to make a distinction among the main lineages.(Suppl. material 1: Fig. S2).

Phylogenetic inference based on both, 3-mtDNA and 15% missing data concatenated datasets, where the same (Fig. 2b; Suppl. material 1: Fig. S1) recovered six main clades with moderate to high levels of statistical support by Maximum-Likelihood and Bayesian inference analyses (Fig. 2b). In both topologies, C. totonacus samples were recovered as a well-supported sister clade of C. ornatus (Fig. 2b) and C. m. molossus was recovered as a sister clade of C. basiliscus, with overall moderate levels of support (Fig. 2b). Also, the clade C. m. nigrescens was recovered forming a well-supported clade apart (Fig. 2b). The most remarkable and novel contribution of our phylogeny is the resolution of the position of the C. m. oaxacus clade, recovered as the basal clade of the group (Fig. 2b; Suppl. material 1: Fig. S1). The network analysis with SplitsTree showed the same pattern, with six well-defined lineages for both dataset (Fig. 2c).

Contrastingly, the haplotype network based only on RAG1 dataset did not show any kind of genealogical structure (Fig. 3): only the C. m. oaxacus lineage yielded two exclusive haplotypes, whereas the remaining haplotypes are randomly distributed over thosputative main clades inferred with both, the 3-mtDNA and 15% missing data concatenated datasets, separately. Most connections between haplotypes contain a considerable number of mutational steps, suggesting the absence of many non-sampled haplotypes. Interestingly, the most abundant, most connected and most central haplotype in the sample is only present in two lineages: C. basiliscus and C. ornatus (Fig. 3).

Genetic distance calculations among six main lineages over the 15% missing data-dataset revealed higher levels of differentiation, ranging from ~5 up to 11% of divergence for all comparisons (Table 1). Considering the average number of nucleotide differences between lineages, the highest value was obtained for C. m. oaxacus vs. C. basiliscus (Dxy= 0.117), whereas the lowest values was inferred for C. totonacus vs. C. ornatus pair (Dxy= 0.075). When the number of net nucleotide substitutions per site between lineages is considered, C. totonacus vs. C. m. oaxacus comparison yielded the highest value (Da= 0.093) and the lowest value of divergence is present in the C. basiliscus vs. C. m. molossus comparison (Da= 0.051).

We selected the three-calibration point-based inference in addition with the 3mtDNA dataset due to its relatively smaller variance in time estimations (Fig. 4a). The divergence between the black-tailed rattlesnakes’ group and his sister group, Crotalus durissus, was dated 8.71 Mya (95% HPD: 5.64–12.62). The split of the most divergent clade within the complex, C. m. oaxacus, was dated 7.71 Mya (95% HPD: 4.99–11.27). The divergence of C. m. nigrescens from the rest of main clades was dated about 7.05 (95% HPD: 4.56–10.41). Divergence of the TMRA of C. totonacus and C. ornatus was dated 5.01 Mya (95% HPD: 2.72–7.79) and the divergence of the TMRA of C. m. molossus and C. basiliscus was dated 6.03 Mya (95% HPD: 3.74–8.85). Reconstruction of the historical demography based on skyride plots indicated a constant population size for condensed lineages molossus-basiliscus (Fig. 4b) and ornatus-totonacus (Fig. 4c) as well as for the C. m. oaxacus clade (Fig. 4e). The reconstruction of C. m. nigrescens indicated remarkable population growth since the last 1.5 Mya (Fig. 4d).

Model selection analysis in BIOGEOBEARS indicates that the best model to reconstruct the ancestral areas of the BEAST tree was DIVALIKE+J (AICc-wt = 0.65). This suggests that the divergence of C. m. oaxacus from the rest of the black-tailed rattlesnakes clades occurred in the Southern Sierra Madre (node: 111, A=45.25%; Fig. 5a, b). Moreover, this area is the most probable ancestral area of the entire group (node: 112, A=15.72%; Fig. 5a, b). The model suggests that the divergence of C. m. nigrescens from the rest of the group occurred, with greater probability, in the TMVB (node: 106, C=24.98%; Fig. 5a, b). The remaining four clades split in the Chihuahuan Desert (node: 83, J=28.41%; Fig. 5a, b). The division between C. ornatus and C. totonacus occurred in this same area (node: 82, J=52.22%; Fig. 5a, b). The ancestral areas of C. ornatus and C. totonacus were the Chihuahuan Desert (node: 79, J=100%; Fig. 5a, b) and the East Sierra Madre (node: 81, K=100%; Fig. 5a, b), respectively, while the area with the highest inferred probability for the MRCA of C. m. molossus and C. basiliscus was the Western Pacific Coastal Plain, Hills, and Canyons (node: 70, G=25.25%; Fig. 5a, b). The ancestral areas of C. m. molossus and C. basiliscus were the Sonoran Desert (node: 69, N=49.88%; Fig. 5a, b) and the Western Pacific Coastal Plain, Hills, and Canyons (node: 65, G=50.91%; Fig. 5a, b), respectively. Dispersal was the predominant biogeographic process explaining the divergence of 46 nodes; the vicariance explained 24 nodes.

Figure 2. 

a. Geographical distribution of the six clades indicated over the topography of North America; b. Phylogenetic inference based on the 3-mtDNA dataset of black-tailed rattlesnakes inferred by Bayesian inference (BI). Nodal support indicates posterior probabilities (BI) above and bootstrap of Maximum likelihood (ML) values below diagonal; c. Neighbor-Net analysis inferred with the SplitsTree method.

Figure 3. 

Minimum Spanning network constructed with RAG1 phased haplotypes. Sizes of colored dots are in agreement with haplotype abundance in the sample and colors represent the main clades recovered along phylogenetic inference (excepting C. totonacus).

Figure 4. 

a. Bayesian divergence time estimation among black-tailed rattlesnake populations based on the 3-mtDNA dataset. Colored clades are indicated as follows: C. basiliscus, green; C. m. molossus, yellow; C. ornatus, red; C. totonacus, purple; C. m. nigrescens, blue; C. m. oaxacus, orange. Skyride reconstructions for the following main condensed lineages; b. C. m. molossus + C. basiliscus; c. C. ornatus + C. totonacus; d. C. m. nigrescens and e. C. m. oaxacus lineages (see the Materials and Methods section).

Figure 5. 

a. Ancestral area reconstruction based on the DIVA+J model from BioGeoBEARS. Pie charts at phylogeny nodes indicate the probability of the ancestral range occupied by that ancestor. Each color represents a reconstructed designated area, the black parts within the pies, indicated the areas with probabilities <5%. The arrows show the vicariance events and the rays the dispersal events. Rhombus (C. basiliscus), circle (C. m. molossus), triangle (C. ornatus), star (C. totonacus), cross (C. m. nigrescens), pentagon (C. m. oaxacus); b. Fifteen North American Terrestrial Ecoregions-Level III considered for the area’s reconstruction; c. Hypothetical scenario explaining the biogeographical structure of the black-tailed rattlesnakes’ group. Photos by Eric Centenero-Alcalá and Sociedad Herpetológica del Noreste de México A.C.

Despite its relevance due its medical importance as well its potential as a biological model organism, several lineages and populations of the rattlesnake genus Crotalus remain poorly understood. In fact, recently, new species have been delimited or even discovered by revisiting long-term known species groups (Bryson et al. 2014; Carbajal-Márquez et al. 2020). Here, we revised the Crotalus molossus complex, considering all the members that comprise the complex and generating information about the historical relationships among their main clades.

Our study is the first to simultaneously evaluate the three recognized members of the Crotalus molossus complex, as well as its previously recognized relatives, C. basiliscus, C. ornatus and C. totonacus, encompassing them in the group of black-tailed rattlesnakes. Both datasets (3-mtDNA and 15% missing data; Fig. 2b, c, Suppl. material 1: Fig. S1) carried enough molecular resolution to support the previous hypothesis provided by Anderson and Greenbaum (2012) and Blair and Sánchez-Ramírez (2016), even though these previous studies included a broader molecular (nuclear) dataset. In fact, they considered C. molossus sensu lato as a paraphyletic species in which C. m. nigrescens was placed in a basal position compared to C. m. molossus-C. basiliscus and C. ornatus-C. totonacus clades (Fig. 1b). One of the main findings in this study was the placement of C. m. oaxacus as the basal clade within the black-tailed rattlesnake group. Only one previous analysis included any sequence of C. m. oaxacus (Wüster et al. 2005). This basal position of C. m. oaxacus was consistently recovered by phylogenetic analyses as well as by network-based approaches. This finding confirms an early divergence of the group in South-Central Mexico (ecoregion A). The latter are in agreement with those suggested by Wüster et al. (2005), who identified South-Central Mexico as the center of origin for the C. durissus group, which has been long recognized as the sister clade of the C. molossus complex.

Our results indicated that the estimated age for the last common ancestor of all black-tailed rattlesnake members dates back to the Miocene at about 7.71 Mya (95% HPD: 4.99–11.27), which is in agreement with the time frame proposed by Anderson and Greenbaum (2012) (7.94 Mya) and by Alencar et al. (2016) (~7.5 Mya), but this was almost twice as old as that proposed by Blair and Sánchez-Ramírez (2016) (~4.8 Mya). The split of C. m. oaxacus was also dated to the Miocene (~7.5) and the divergence of C. m. molossus and C. basiliscus occurred during the Pliocene, while the most recent divergence events dated within the group, calculated for C. m. nigrescens, occurred approximately <1 Mya (Fig. 4a).

Historical demographic reconstruction indicated that most lineages yielded “constant” historical population sizes (Fig. 4b, c, e), suggesting that population expansions were not typical along the divergence of this group, at least for mitochondrial markers. We were only able to find evidence of population growth in C. m. nigrescens, occurred within a short and relatively recent time frame (~1.5 Mya) (Fig. 4d). Considering that this clade is the most ecologically widespread within this group, this recent population growth could be consistent with an expansion scenario (see below). A wider sampling from the rest of the members is required with the aim of ensuring that this expansion pattern is exclusively for this lineage.

While Anderson and Greenbaum (2012) previously suggested a basal divergence of the C. molossus group, explained by a Neogene vicariance event at the SMOcc, our results placed the origin and early divergence of the C. molossus complex in the Sierra Madre de Sur (ecoregion A), splitting ~9 Mya from its sister lineage, the Neotropical rattlesnake group C. durissus. We infer that the ancestral area of this group was the northwestern part of the Sierra Madre de Sur (ecoregion A), in contact with the TMVB (ecoregion C). This divergence was likely related with the early orogenetic process of the TMVB (Fig. 5a–c). Deep divergences in other co-distributed vertebrates have been broadly attributable to the eastern uplift of the TMVB, splitting lineages between Central Mexico and the Oaxaca highlands (León-Paniagua and Morrone 2009; Bryson et al. 2011a, b; Bryson et al. 2014). Our reconstruction indicated that the northeastern portion of the Sierra Madre de Sur (ecoregion A) as well as the eastern TMVB (ecoregion C) could be plausible regions of early divergence and northward dispersal, both starting ~7.9 Mya (Fig. 5b, c). Black-tailed rattlesnakes are generalists, largely associated with mountains and foothills, but specimens have been recorded from sand dunes and xeric scrubs, suggesting highly dispersing capabilities (Anderson and Greenbaum 2012).

The subsequent split, which comprises divergence of C. m. nigrescens was inferred to happen at the northwest area of Jalisco (ecoregions C, E and H, Fig. 5a, b). Our reconstruction indicated a population expansion towards the east, in agreement with our Skyride reconstruction, which promoted its divergence into three main subclades: one distributed across SMOcc pine and oak forest (ecoregion H), another occupying the TVBM (ecoregion C) and another distributed across Lerma depression (E, Fig. 5b, c). Several phylogeographical studies along the frontier between Lerma depression and TMVB have suggested the effect of Pliocene and Pleistocene geological dynamics promoting multiple events of divergence (Huidrobo et al. 2006; Domínguez-Domínguez et al. 2008; Hernández-Leal et al. 2019).

Northward to the TMVB, our estimations indicated that the split between the ancestors of C. m. molossus-C. basiliscus and C. ornatus-C. totonacus clades, occurred likely in the Chihuahuan Desert (Fig. 5a–c). The latter agrees with that described by Anderson and Greenbaum (2012); however, in this region there is a narrow desert corridor that allows for exchange of taxa between these deserts, named the Cochise Barrier Filter, where other authors have found taxa with distributions on both sides of the Cochise Barrier Filter (Pyron and Burbrink 2010; Myers et al. 2017; Myers et al. 2019). In this region, the Mimbres River and the Western Continental Divide, which is north of the Mimbres Valley, seem to represent the only potential barrier that may explain the apparent lack of introgression between C. ornatus and C. molossus. To explain the split between the C. basiliscus and C. m. molossus clades, we hypothesized that its ancestor, distributed in the southern Chihuhuan Desert (ecoregion J), colonized the coastal plains of the Pacific Coast (ecoregion G), intersecting the SMOcc (ecoregion H) and arriving in the Western Pacific Plain and Hills (ecoregion F). These corridors probably correspond to the northwestern area of Jalisco, where some authors have indicated the occurrence of a secondary contact zone between C. basiliscus and C. m. nigrescens (Klauber 1952; McCranie 1981).

From this central area of the Pacific Coast (F), some populations dispersed to the south (ecoregion D) and north (ecoregion G) areas, where those ancestral populations gave rise to the C. basiliscus clade. This clade could have occupied tropical deciduous forests, naturally present along the Mexican Pacific Coast, as a corridor to spread north and south. Certain populations that headed north toward the Sonoran Desert (ecoregion N) later originated in the Crotalus m. molossus clade, which then spread across the ecoregions of the Madrean Archipelago (ecoregion M) and the Arizona/New Mexico Mountains (ecoregion O). The northernmost populations of C. basiliscus probably differed from the populations that gave rise to C. m. molossus due to the local adaptation to the different habitats (Fig. 5)

Alternatively, the divergence of C. molossus- C. bassiliscus ancestors may have occurred further north after its ancestral populations reached their northernmost extent at the SMOcc in Sonora Desert (ecoregion N). This split may had been ecologically driven by vegetation replacement, from the Oligocene to the late Miocene, when subtropical thorn forests and woodland savannas were gradually replaced with savanna and semidesert habitats, and this continued during the latest Miocene with the expansion of regional deserts, grasslands and shrub-steppes (Alexander and Riddle 2005). Subsequently, C. basilicus spread along the Tropical deciduous forest, present widely along the Mexican Pacific (Fig. 5a–c); this divergence pattern has been reported in snakes distributed along the Mexican Pacific Coast (Trimorphodon, Devitt 2006).

Finally, a similar pattern may explain the divergence and distribution of C. totonacus with respect to C. ornatus; the Totonacan Rattlesnake is endemic to Mexico, where it inhabits the coastal plain of the Gulf of Mexico and adjacent areas on the eastern versant of the SMOr, and has been reported from the states of Nuevo León, Tamaulipas, San Luis Potosí, Querétaro, Hidalgo, and Veracruz (Ramírez-Bautista et al. 2014; Farr et al. 2015; Ramos-Frías et al. 2015). Divergence between the ancestral populations of those lineages would take place in response to continuing isolation of the regional deserts throughout the late Miocene, early Pliocene and more recently, during glacial cycles of the Pleistocene (Alexander and Riddle 2005; Bryson et al. 2011c).

Wüster et al. (2005) demonstrated that C. molossus was highly heterogeneous, suggesting that it may also represent a species complex. Anderson and Greenbaum (2012) noted unexpected diversity within this group and inferred that the C. m. nigrescens lineage is the basal member of the group. However, in the present study, we included C. m. oaxacus sequences and established that this is the basal clade within the group. Anderson and Greenbaum (2012) also determined a divergence between the USA western lineage C. m. molossus and the eastern formally erected species C. ornatus (Hallowell, 1854). In our study, we confirmed that the divergence between C. m. molossus and C. basiliscus and between C. ornatus and C. totonacus (Alencar et al. 2016; Blair and Sánchez-Ramírez 2016). While C. ornatus, C. basiliscus and C. totonacus have been traditionally considered species, our genealogical inference provides evidence that the three subspecies within the complex could be also considered species under a phylogenetic concept (De Queiroz 2007).

Our findings suggest that three lineages within the Crotalus molossus complex, which are in agreement with previously considered subspecies, may be potentially recognized as species, like in the case of Crotulus ornatus and Crotalus basiliscus tacking into account their significant levels of genetic divergence (5–11%) as well as their remarkable geographical structure among them: C. m. molossus (Northern black-tailed rattlesnake) (Baird and Girard 1853) distributed from Sonora, Chihuahua and Coahuila to Central Texas, New Mexico, and Arizona; C. m. nigrescens (Mexican black-tailed rattlesnake) (Gloyd, 1936) distributed through Central Mexico at north of TMVB as far north as west-central Chihuahua and C. m. oaxacus (Oaxacan black-tailed rattlesnake) (Gloyd, 1948) in Oaxaca and southeastern Puebla.

On the other hand, Ruíz-Sánchez et al. (2019) found that Crotalus estebanensis (Klauber, 1949) was the sister to the C. basiliscus from Nayarit, Mexico (based on mitochondrial regions 12S and 16S), with an estimated divergence of ~3.22 Ma, which is consistent with the origin of San Esteban Island during the Miocene and the Pliocene. Future analyses considering wider samples of C. basiliscus and including common genes for all main clades, will provide enough resolution to clarify if the C. estebanensis population represents a completely differentiated clade from C. basiliscus or if C. estebanensis haplotypes are nested within the C. basiliscus lineage.

We are aware that analyses based solely on mtDNA provide a restricted point of view of species evolutionary history (i.e., matrilineal inheritance) and unfortunately, availability of nuclear data from vipers is still scarce. While maternally inherited information has been long recognized to lead to significant information, such as biogeographical discoveries (Riddle et al. 2000; McGuire et al. 2007) and the occurrence of cryptic species (Zink and Barrowclough 2008; Zink 2010), all these results should be taken with caution as our analysis is based mainly on mitochondrial genes, and hybridization or mitochondrial introgression could be responsible for those patterns (see Reyes-Velasco et al. 2022). Future analyses considering a wider nuclear loci sampling can test whether mtDNA-based genealogy accurately reflects the species phylogenic relationships and it will provide insights into other processes, such as hybridization among C. molossus members.