Research Article
Research Article
A tale about vipers’ tails: phylogeography of black-tailed rattlesnakes
expand article infoVíctor Hugo Muñoz-Mora, Marco Suárez-Atilano§, Ferruccio Maltagliati|, Fabiola Ramírez-Corona, Alejandro Carbajal-Saucedo#, Ruth Percino-Daniel¤, Joachim Langeneck|, Maristella D’Addario«, Armando Sunny«
‡ University of Ferrara, Ferrara, Italy
§ University of California, Berkeley, United States of America
| Università di Pisa, Pisa, Italy
¶ Universidad Nacional Autónoma de México, Ciudad de México, Mexico
# ISBI Biosciences S.C., Huitzilac, Mexico
¤ Universidad de Ciencias y Artes de Chiapas, Tuxtla Gutiérrez, Mexico
« Universidad Autónoma del Estado de México, Toluca, Mexico
Open Access


The phylogenetic relationships among black-tailed rattlesnakes remain poorly understood and some authors indicated that the diversity of this group has been underestimated and additional analyses are required to clarify the biogeographic patterns throughout its distribution in Mexico. Therefore, the aim of this study was to elucidate the phylogenetic relationships among black-tailed rattlesnakes across their range, identifying relative divergence times among the main clades and reconstructing the biogeographical history of the group. Three partial mitochondrial genes (ND4, cytb and ATPase6) and one nuclear gene (RAG1) were sequenced to infer the phylogenetic relationships, through the maximum likelihood and Bayesian inference-based methods; demographic history reconstruction was investigated through Bayesian Skyline plot analysis and the ancestral area reconstruction was carried out considering a Bayesian framework. We found strong evidence that the black-tailed rattlesnakes’ group is composed of six clades, which is in agreement with subspecies previously reported. Divergence time estimation indicated that the origin of the C. molossus group could be traced to the middle of the Miocene (~7.71 Mya). Ancestral area reconstruction indicated that early divergence events occurred in Central Mexico, probably related to the geological dynamics of the Trans-Mexican Volcanic Belt. The lineage C. m. oaxacus is the basal member of the C. molossus group. Furthermore, the combination of geological events and changes in Quaternary vegetation may have contributed to the divergence of C. molossus clades. Our results suggest several clades within C. molossus complex could be potentially recognized as separate species.

Key Words

ancestral area reconstruction, Crotalus molossus, divergence time estimation, Mexico, molecular biogeography, North America, phylogenetics, phylogeography, Viperidae


Phylogenetics is fundamental in evolutionary biology to elucidate diversification over time and space, and it is irreplaceable in understanding the evolutionary history of full biotas (Riddle et al. 2000; Graham et al. 2018). Phylogenetic hypotheses have been used to infer associations between the timing of lineage divergence and past geologic and climatic events from Earth’s history (Mantooth and Riddle 2011), through approaches such as molecular biogeography and phylogeography (Riddle et al. 2008; Mantooth and Riddle 2011; Graham et al. 2018). The application of phylogenetic analyses is not restricted only to clarify the evolutionary relationships between species or higher taxa, when populations or individuals within species are considered as Operational Taxonomic Units (OTUs); different evolutionary processes can be detected and evaluated such as gene flow, hybridization and demographic changes among populations (Avise 2000). If we also considered a biogeographical dimension in this kind of inference, we would be able to elucidate how the distribution of species or populations has changed over the time (Freeland 2005).

The physiography of Mexico has been shaped by a dynamic combination of geological, climatic, and historical processes, and its resulting complex topography has provided a matrix for the evolution of a spectacularly diverse biota characterized by high endemism (Bryson et al. 2011). Multiple biotic assemblages are likely to have been shaped along the frontier between the Neotropics and the Nearctic, which is the Mexican Transition Zone (MTZ) located in Central Mexico, and they are associated with major geographical features, such as the Sierra Madre Occidental (SMOcc), the Trans-Mexican Volcanic Belt (TMVB), and the Sierra Madre Oriental (SMOr) (Savage 1982).

Rattlesnakes (Crotalus and Sistrurus) constitute a diverse group of endemic New World vipers, ranging from Canada to Argentina (Klauber 1972) and occupying a wide spectrum of ecosystems, including deserts, cloud forests and alpine regions (Campbell and Lamar 2004). Rattlesnakes had received considerable attention, mainly due to their medical importance and broad geographical and ecological distribution (Blair and Sánchez-Ramírez 2016). In addition, these vipers are considered key in the control of rodent populations, making them a priority group for conservation (Valencia-Hernández 2006). Mexico has the highest diversity of rattlesnakes (Campbell and Lamar 2004; Place and Abramson 2004; Uetz et al. 2021), including 43 of the 53 described species in the genus Crotalus, and 26 of them are endemic (Uetz et al. 2021).

The central portion of the country has been considered an important diversification region for this genus, considering the presence of many basal groups in the area, such as the high-mountain rattlesnakes C. polystictus (Cope, 1865) and C. ravus (Cope, 1865), as well the C. triseriatus (Wagler, 1830) complex (Murphy et al. 2002; Campbell and Lamar 2004; Bryson et al. 2011). Previous approaches have indicated the pine-oak forest of Central Mexico as the most likely ancestral habitat for the genus Crotalus (Murphy et al. 2002; Campbell and Lamar 2004; Place and Abramson 2004; Blair and Sánchez-Ramírez 2012). Within this genus, Crotalus molossus (Baird & Girard, 1853) has one of the widest ecological and geographical distributions (Fig. 1), ranging from western, central, and southern Arizona, central and southern New Mexico, and southwestern and central Texas in the United States, south through Mexico to the southern edge of the Mexican Plateau and Mesa del Sur (Oaxaca), occupying regions from Sonora, Chihuahua, and Coahuila through the Mesa Central to the Southern Sierra Madre in Oaxaca (Fig. 1a) and inhabiting a range of biomes, from pine-oak forests to xerophilous scrub, at elevations from sea level to 2930 m a.s.l. (Stebbins 2003; Campbell and Lamar 2004; Ramírez-Bautista et al. 2009).

Figure 1. 

a. Geographical distributions of C. m. molossus (yellow), C. m. nigrescens (blue), C. m. oaxacus (orange) as well as the long-recognized species C. basiliscus (green), C. ornatus (red), and C totonacus (purple) [modified from McCranie (1981) and Anderson and Greenbaum (2012)]; b. Phylogenetic relationships among black-tailed rattlesnakes, according to the hypothesis proposed by Wüster et al. (2005) and Anderson and Greenbaum (2012). While relationships between northern C. molossus members have been analyzed, the phylogenetic position of C. m. oaxacus within the complex is unknown.

Traditionally, three subspecies have been recognized: Crotalus molossus molossus (Baird and Girard 1853), Crotalus molossus nigrescens (Gloyd 1936), and Crotalus molossus oaxacus (Gloyd 1948; Crotalus molossus sensu lato; Anderson and Greenbaum 2012). Controversially, phylogenetic relationships within this group remain poorly understood (Anderson and Greenbaum 2012; Blair and Sánchez-Ramírez 2016). C. molossus was initially placed within Crotalus durissus (Linnaeus, 1758) group (Gloyd 1940; Klauber 1956; Brattstrom 1964; Foote and MacMahon 1977), later it was identified as a sister clade of C. atrox (Baird and Girard 1853). Recent studies with C. basiliscus (Linnaeus, 1758) support a closer relationship (Murphy et al. 2002; Wüster et al. 2005; Castoe and Parkinson 2006; Blair and Ramírez-Sánchez 2016).

The phylogeographic structure of the northern populations of the C. molossus group was analyzed by Anderson and Greenbaum (2012), who found deep genetic and morphological differentiation between southern USA populations, concluding that eastern populations are a separate species: Crotalus ornatus (Hallowell 1854). Furthermore, they found support for sister relationships between Crotalus totonacus (Gloyd and Kauffeld 1940) and the eastern C. ornatus lineage, as well as between C. basiliscus and the western C. m. molossus lineage (Fig. 1b). However, the same authors also indicated that the diversity of this group has been underestimated and additional analyses are required to clarify the biogeographic patterns throughout its distribution in Mexico, since the phylogenetic relationships between C. basiliscus and the southernmost subspecies, C. m. nigrescens and C. m. oaxacus, are still poorly understood (Fig. 1b).

Therefore, the aims of this study were the following: a) to clarify the phylogenetic relationships within the black-tailed rattlesnake group, b) to estimate the temporary framework in which the main divergence events among clades occurred, and, finally, c) to reconstruct the most likely ancestral distribution of those taxa to identify the main biogeographical events responsible for their present divergence and distribution. In order to clarify the taxonomical perspectives in this study, we defined the following six main lineages as part of the black-tailed rattlesnake group: C. m. molossus, C. m. nigrescens, C. m. oaxacus, C. ornatus, C. basiliscus, and C. totonacus.

Materials and methods

Sampling and molecular procedures

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 MN327706MN327962, 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 ddH2O, 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 (, to ensure correct identification of amplifications. Subsequently, the sequences were traduced to proteins to inspect the presence of unexpected stop codons.

Phylogenetic inference

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.

Genetic diversity and population structure analyses

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.

Divergence time estimation and historical demographic reconstruction

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 1000th 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 107 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 1010 generations, discarding 10% as burn-in and with coalescent intervals (m) set to 5.

Ancestral area reconstruction

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).


Phylogenetic inference and genetic diversity

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 FS tests (FS= –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).

Divergence time estimation and historical demographic reconstruction

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).

Ancestral areas reconstruction

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).

Table 1.

Genetic distances matrix estimated considering the 15% missing data-dataset, among six main clades. The average number of nucleotide substitutions per site between populations (Dxy) values are below diagonal and the number of net nucleotide substitutions per site between populations (Da) above diagonal. In all cases standard deviation values are within parentheses.

Dxy /Da Crotalus basiliscus C. m. molossus C. ornatus C. totonacus C. m. niegrescens C. m. oaxacus
Crotalus basiliscus 0.05079 (0.02526) 0.05406 (0.01975) 0.07017 (0.03333) 0.07006 (0.0152) 0.07595 (0.02928)
C. m. molossus 0.10132 (0.02492) 0.05495 (0.01500) 0.07417 (0.02818) 0.07523 (0.01194) 0.08231 (0.02565)
C. ornatus 0.10159 (0.01951) 0.07901 (0.01463) 0.06330 (0.01743) 0.06899 (0.00800) 0.07519 (0.01688)
C. totonacus 0.10902 (0.03329) 0.08954 (0.02805) 0.07568 (0.01731) 0.09229 (0.01536) 0.09367 (0.03140)
C. m. niegrescens 0.11241 (0.0149) 0.09411 (0.01167) 0.08487 (0.00769) 0.09949 (0.01534) 0.09060 (0.01255)
C. m. oaxacus 0.11773 (0.02919) 0.10061 (0.02551) 0.09049 (0.01671) 0.10029 (0.03138) 0.10072 (0.01248)
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.


Phylogenetics analyses

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.

Divergence time estimations and ancestral area reconstruction

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).

Taxonomical implications

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.


We are deeply grateful to Dr. Wolfgang Wüster for his valuable comments and suggestions on improving the manuscript. V.H.M.M. is grateful to the graduate program Maestría en Ciencias Agropecuarias y Recursos Naturales at Universidad Autónoma del Estado de México (UAEM: 0711983) and Consejo Nacional de Ciencia y Tecnología scholarship (CONACyT: 785320). Likewise, we thank all those who helped in the tissue collection, especially Agustín Álvarez-Trillo. We extend our appreciation to Emma Taylor for language reviewing. We also thank Eric Centenero-Alcalá and Sociedad Herpetológica del Noreste de México A.C. for kindly sharing his images of rattlesnakes. This study was conducted with the approval of the ethics committee of Universidad Autónoma del Estado de México and the field permit: SGPA/DGVS/011587/17 SEMARNAT. A.S received financial support from Secretaría de Investigación y Estudios Avanzados (SYEA) of Universidad Autónoma del Estado de México (Grant: 4732/2019CIB). This paper was completed while M.S.A was on his Postdoctoral stay at the Museum of Vertebrate Zoology at UC Berkeley (UCMEXUS-CONACyT, 10327268). The authors declare that they have no conflict of interest.


  • Alencar LRV, Quental TB, Grazziotin FG, Alfaro M, Martins M, Venzon M, Zaher H (2016) Diversification in vipers: Phylogenetic relationships, time of divergence and shifts in speciation rates. Molecular Phylogenetics and Evolution 105: 50–62.
  • Anderson CG, Greenbaum E (2012) Phylogeography of northern populations of the black-tailed rattlesnake Crotalus molossus Baird & Girard, 1853, with the revalidation of C. ornatus Hallowell, 1854. Herpetological Monographs 26: 19–57.
  • Baird SF, Girard C (1853) Catalogue of North American Reptiles in the Museum of the Smithsonian Institution. Part 1.-Serpents. Smithsonian Inst., Washington, 172 pp.
  • Blair C, Sánchez-Ramírez S (2016) Diversity-dependent cladogenesis throughout western Mexico: Evolutionary biogeography of rattlesnakes Viperidae: Crotalinae: Crotalus and Sistrurus. Molecular Phylogenetics and Evolution 97: 145–154.
  • Bryson Jr RW, Murphy RW, Lathrop AL, Lazcano-Villareal DL (2011a) Evolutionary drivers of phylogeographical diversity in the highlands of Mexico: A case study of the Crotalus triseriatus species group of montane rattlesnakes. Journal of Biogeography 38: 697–710.
  • Bryson Jr RW, García-Vázquez U, Riddle BR (2011b) Phylogeography of Middle American gophersnakes: mixed responses to biogeographical barriers across the Mexican Transition Zone. Journal of Biogeography 38: 1570–1584.
  • Bryson Jr RW, Murphy RW, Graham MR, Lathrop A, Lazacano D (2011c) Ephemeral Pleistocene woodlands connect the dots for highland rattlesnakes of the Crotalus intermedius group. Journal of Biogeography 38: 2299–2310.
  • Bryson Jr RW, Linkem CW, Dorcas ME, Lathrop A, Jones JM, Alvarado-Díaz J, Grünwald CI, Murphy RW (2014) Multilocus species delimitation in the Crotalus triseriatus species group Serpentes: Viperidae: Crotalinae, with the description of two new species. Zootaxa 38233: 475–496.
  • Campbell JA, Lamar WW (2004) Venomous reptiles of the Western Hemisphere. Cornell University Press, Ithaca, New York, 528 pp.
  • Cope ED (1885) Twelfth contribution to the herpetology of tropical America. Proceedings of the American Philosophical Society 22: 167–194.
  • Darriba D, Taboada G, Doallo R, Posada D (2012) jModelTest 2: more models, new heuristics and parallel computing. Nature Methods 9: 772.
  • Devitt TJ (2006) Phylogeography of the Western Lyresnake Trimorphodon biscutatus: testing arid land biogeographical hypotheses across the Nearctic–Neotropical transition. Molecular Ecology 15: 4387–4407.
  • Domínguez-Domínguez O, Alda F, De León GPP, García‐Garitagoitia J, Doadrio I (2008) Evolutionary history of the endangered fish Zoogoneticus quitzeoensis Bean, 1898 Cyprinodontiformes: Goodeidae using a sequential approach to phylogeography based on mitochondrial and nuclear DNA data. BMC Evolutionary Biology 81: 161.
  • Douglas ME, Douglas MR, Schuett GW, Porras LW (2006) Evolution of rattlesnakes Viperidae; Crotalus in the warm deserts of western North America shaped by Neogene vicariance and Quaternary climate change. Molecular Ecology 15: 3353–3374.
  • Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology and Evolution 29: 1969–1973.
  • Earl D, VonHoldt BM (2014) Structure harvester. Version 0.6, 94: 2014.
  • Farr WL, Nevárez De los Reyes M, Banda-Leal J, Lazcano D (2015) The distribution of Crotalus totonacus in Nuevo León, Mexico. Mesoamerican Herpetology 2: 243–251.
  • Gloyd HK (1936) A Mexican subspecies of Crotalus molossus Baird and Girard. Occasional Papers of the Museum of Zoology, University of Michigan 325: 1–5.
  • Gloyd HK (1940) The Rattlesnakes, Genera Sistrurus and Crotalus. A Study in Zoogeography and Evolution. Chicago, Chicago Academy of Sciences, 270 pp.
  • Gloyd HK (1948) Description of a neglected subspecies of rattlesnake from Mexico. Natural History Miscellanea 17: 1–4
  • Gloyd HK, Kauffeld CF (1940) A new rattlesnake from Mexico. Bulletin of the Chicago Academy of Sciences 6(2): 11–14.
  • Guiher TJ, Burbrink FT (2008) Demographic and phylogeographic histories of two venomous North American snakes of the genus Agkistrodon. Molecular Phylogenetics and Evolution 482: 543–553.
  • Guindon S, Dufayard JF, Lefort V, Asimova M, Hordijk W, Gascuel O (2010) New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systematic Biology 59: 307–321.
  • Hallowell E (1854) Descriptions of new reptiles from California. Proceedings of the Academy of Natural Sciences, Philadelphia 7: 91–97.
  • Hernández-Leal MS, Suárez-Atilano M, Piñero D, González-Rodríguez A (2019) Regional patters of tenetic structure and environmetal differentiation in willow populations Salix humboldtiana Wild. from Central Mexico. Ecology and Evolution 9: 9564–9579.
  • Huidobro L, Morrone JJ, Villalobos JL, Álvarez F (2006) Distributional patterns of freshwater taxa fishes, crustaceans and plants from the Mexican Transition Zone. Journal of Biogeography 334: 731–741.
  • Klauber LM (1972) Rattlesnakes: their habitats, life histories, and influences on mankind, 2 volumens. 2nd Edn. University of California Press, Berkeley, 400 pp.
  • Landis MJ, Matzke NJ, Moore BR, Huelsenbeck JP (2013) Bayesian Analysis of Biogeography when the Number of Areas is Large, Systematic Biology 626: 789–804.
  • León-Paniagua L, Navarro-Sigüenza AG, Hernández-Baños BE, Morales JC (2007) Diversification of the arboreal mice of the genus Habromys Rodentia: Cricetidae: Neotominae in the Mesoamerican highlands. Molecular Phylogenetics and Evolution 42: 653–664.
  • León-Paniagua LL, Morrone JJ (2009) Do the Oaxacan Highlands represent a natural biotic unit? A cladistics biogeographical test based on vertebrate taxa. Journal of Biogeography 36: 1939–1944.
  • Linnaeus C (1758) Systema naturæ per regna tria naturæ, secundum classes, ordines, genera, species, cum characteribus, differentiis, synonymis, locis. Tomus I. Editio decima, reformata. Laurentii Salvii, Holmiæ. 10th Edn., 824 pp.
  • Linnaeus C (1766) Systema naturæ per regna tria naturæ, secundum classes, ordines, genera, species, cum characteribus, differentiis, synonymis, locis. Tomus I. Editio duodecima, reformata. Laurentii Salvii, Stockholm, Holmiae, 1–532.
  • Matzke-Nicholas J (2013) Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Frontiers of Biogeography 54: 242–248.
  • McCraine JR (1981) Crotalus basiliscus Cope. Mexican west coast rattlesnake. Catalogue of American Amphibians and Reptiles 283: 1–2.
  • McGuire JA, Linkem CW, Koo M, Hutchison MDW, Lappin AK, Orange DO, Lemos-Espinal J, Riddle BR, Jaeger J (2007) Mitochondrial introgression and incomplete lineage sorting through space and time: Phylogenetics of crotaphytid lizards. Evolution 6112: 2879–2897.
  • Murphy RW, Fu J, Lathrop A, Feltham JV, Kovak V (2002) Phylogeny of the rattlesnakes Crotalus and Sistrurus inferred from sequences of five mitochondrial DNA genes. Biology of the vipers, 69–92.
  • Myers EA, Hickerson MJ, Burbrink FT (2017) Asynchronous diversification of snakes in the North American warm deserts. Journal of Biogeography 442: 461–474.
  • Myers EA, Xue AT, Gehara M, Cox C, Davis-Rabosky AR, Lemos-Espinal J, Martínez-Gómez JE, Burbrink FT (2019) Environmental Heterogeneity and Not Vicariant Biogeographic Barriers Generate Community Wide Population Structure in Desert Adapted Snakes. Molecular Ecology 2820: 4535–4548.
  • Navarro-Sigüenza AG, Peterson AT, Nyari A, García-Deras G, García-Moreno J (2008) Phylogeography of the Buarremon brush-finch complex Aves, Emberizidae in Mesoamerica. Molecular Phylogenetics and Evolution 47: 21–35.
  • Parkinson CL, Campbell JA, Chippindale PT, Schuett G (2002) Multigene phylogenetic analysis of pitvipers, with comments on their biogeography. Biology of the Vipers 9: 3–110.
  • Parkinson CL, Zamudio KR, Greene HW (2000) Phylogeography of the pitviper clade Agkistrodon: historical ecology, species status, and conservation of cantils. Molecular Ecology 94: 411–420.
  • Place AJ, Abramson CI (2004) A Quantitative Analysis of the Ancestral Area of Rattlesnakes. Journal of Herpetology 381: 151–156.
  • Pyron RA, Burbrink FT (2010) How does ecological opportunity influence rates of speciation, extinction, and morphological diversification in New World rattlesnakes tribe Lampropeltini?. Evolution 644: 934–943.
  • Rafinesque CS (1818) Further accounts of discoveries in natural history in the western states. American Monthly Magazine and Critical Review 4(1): 41.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Systematic Biology, syy032.
  • Ramírez-Bautista A, Hernández-Salinas U, García-Vázquez UO, Leyte-Manrique A, Canseco-Márquez L (2009) Herpetofauna del Valle de México: Diversidad y Conservación. Universidad Autónoma del Estado de Hidalgo, Hidalgo, México.
  • Ramos-Frías J, Fernández-Badillo L, Aguilar-López M, Ángeles-Escudero JI (2015) Geographic Distribution. Crotalus totonacus Totonacan Rattlesnake. Herpetological Review 46: 219.
  • Ree RH, Smith SA (2008) Maximum likelihood inference of geographic range evolution by dispersal, local extinction and cladogenesis. Systematic Biology 57: 4–14.
  • Riddle BR, Dawson MN, Hadly EA, Hafner DD, Hickerson MJ, Mantooth SJ, Yoder AD (2008) The role of molecular genetics in sculpting the future of integrative biogeography. Progress in Physical Geography 32: 173–202.
  • Riddle BR, Hafner DJ, Alexander LF, Jaeger JR (2000) Cryptic vicariance in the historical assembly of a Baja California Peninsular Desert biota. Proceedings of the National Academy of Sciences 9726: 14438–14443.
  • Ruíz-Sánchez E, Arnaud G, Cruz-Andrés OR, García-De León FJ (2019) Phylogenetic relationships and origin of the rattlesnakes of the Gulf of California islands Viperidae: Crolalinae: Crotalus. Herpetological Journal 29: 162–172.
  • Savage JM (1982) The enigma of the Central American herpetofauna: dispersal or vicariance? Annals of the Missouri Botanical Garden 69: 464–547.
  • Stebbins RC (2003) A Field Guide to Western Reptiles and Amphibians. Third Edition. Houghton Mifflin Company, Boston, Massachusetts.
  • Stephens M, Donnelly P (2003) A comparison of Bayesian methods for haplotype reconstruction from population genotype data. American Journal of Human Genetics 73: 1162–1169.
  • Valencia-Hernández AA (2006) Taxonomía y distribución del género Crotalus Linneo, 1758 en el Estado de Hidalgo. Bachelor’s degree Thesis, Universidad Autónoma del Estado de Hidalgo, Hidalgo, México, 114 pp.
  • Velo-Antón G, Parra J. Parra-Olea G, Zamudio KR (2013) Tracking climate change in a dispersal-limited species: reduced spatial and genetic connectivity in a montane salamander. Molecular Ecology 22: 3261–3278.
  • Wagler JG (1830) Natürliches System der Amphibien, mit vorangehender Classification der Säugetiere und Vögel. Ein Beitrag zur vergleichenden Zoologie. 1.0. Cotta, München, Stuttgart, and Tübingen, 354 pp.
  • Wilken E, Jiménez-Nava F, Griffith G (2011) North American Terrestrial Ecoregions-Level III. Commission for Environmental Cooperation, Montreal, Canada, 149 pp.
  • Wüster W, Ferguson JE, Quijana-Mascareñas JA, Pook CE, Salomão MG, Thorpe RS (2005) Tracing an invasion: landbridges, refugia and the phylogeography of the Neotropical rattlesnake Serpentes: Viperidae: Crotalus durissus. Molecular Ecology 14: 1095–1108.
  • Yu Y, Blair C, He X (2020) RASP 4: ancestral state reconstruction tool for multiple genes and characters. Molecular Biology and Evolution 37(2): 604–606.

Supplementary material

Supplementary material 1 

Tables S1–S4, Figures S1–S5

Víctor Hugo Muñoz-Mora, Marco Suárez-Atilano, Ferruccio Maltagliati, Fabiola Ramírez-Corona, Alejandro Carbajal-Saucedo, Ruth Percino-Daniel, Joachim Langeneck, Maristella D'Addario, Armando Sunny

Data type: Occurences and Genbank ID, phylogenetic analysis and images

Explanation note: Occurrences and Genbank ids, as well as phylogenetic analyzes carried out that complement the study.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (1.75 MB)
login to comment