In total we measured 684 specimens (Table 1) deposited at the Colección Nacional de Anfibios y Reptiles (CNAR), ENCB and Museo de Zoología Alfonso L. Herrera (MZFC). Additionally, to verify some historical records and for comparison and objective reference, we requested photographs of museum specimens including type material deposited in another 19 collections (See supplementary file 1: Museum specimens).

We georeferenced all localities using GoogleEarth Pro v.7.3.3.7699 and digitized topographic maps available in the digital library of the Instituto Nacional de Estadística y Geografía (INEGI, https://www.inegi.org.mx/app/mapas). In the field we used a Garmin etrex30 GPS with WGS84 datum to record collection localities.

For genetic analyses, we obtained 56 tissue samples from the MZFC collection and field work, that include individuals collected in close proximity to the type localities of all five recognized subspecies of S. torquatus, as well as the undescribed Sceloporus sp. from western Mexico sensu Martínez-Méndez and Méndez-De la Cruz (2007). Samples of Sceloporus bulleri Boulenger, 1895, S. mucronatus Cope, 1885, and Sceloporus grammicus Wiegmann, 1828 were also included (Fig. 1; Supplementary file 2: Tissue sampling and GenBank accession numbers). We chose the mitochondrial 12S and ND4 loci, and the nuclear RAG1 locus for genetic analyses, as these regions have successfully been utilized to delimit species Sceloporus species in similar studies (Wiens and Penkrot 2002; Martínez-Méndez et al. 2012; Díaz-Cárdenas et al. 2017).

To perform DNA extractions, we used the Qiagen™ DNeasy Blood & Tissue Kit™ following the manufacturer’s protocol.

We amplified fragments of the 12S and ND4 mtDNA regions, and RAG1 of nDNA by means of polymerase chain reaction (PCR) under the following standardized conditions: 1μL DNA extraction, 9.45μL dH₂O, 3μL 5X MyTaq™ Reaction Buffer, 0.5µL Primer F [10µM], 0.5µL Primer R (10µM) and 0.15µL MyTaq™ Bioline™ (5U). PCRs were carried out in a Multigene Optimax LabNet™ thermocycler with the following annealing temperatures for each molecular marker: 45°C, 12S; 54°C, ND4; and 50°C, RAG1. The oligonucleotides sequences used (Table 2) were taken from Kocher et al. (1989), Forstner et al. (1995) and Wiens et al. (2010).

We used the sequencing service of the Laboratorio Nacional de Biodiversidad (LANABIO) at the Instituto de Biología (IBUNAM), which uses the BigDye Terminator v.3.1 Applied Biosystems kit and a final purification with Sephadex G-50 before analyzing cycle sequencing product on an Applied Biosystems 3730 xL DNA Analyzer Sequencer.

Once sequences were obtained, we used MUSCLE (Edgar 2004) implemented in MEGA-X v.10.0.5 (Kumar et al. 2018) to pair contigs and align sequences. Subsequently, we reviewed alignments by eye, and eliminated small regions of the sequences that contained polymorphic sites that were difficult to align.

We constructed two molecular data matrices —the first one exclusively with the mtDNA data (12S + ND4) and the second with the combined data from mtDNA + nDNA (12S + ND4 + RAG1). We also included sequence data generated in previous works (Martínez-Méndez and Méndez-De la Cruz 2007; Leaché and Mulcahy 2007). For accession numbers of sequences used see supplementary file 2: Tissues sampling. To identify the optimal partitions in both datasets, as well as the best nucleotide substitution model for each partition, we used PartitionFinder2 (Lanfear et al. 2016) through the CIPRES Science Gateway v.3.3 interface (Miller et al. 2010), with potential partitions divided by codon position for coding regions.

To infer the phylogenetic relationships of S. torquatus ssp. we performed both Bayesian inference and ML analyses with both mitochondrial and combined datasets, using MrBayes v.3.2.7a (Ronquist et al. 2012) and RaxML-HPC2 (Stamatakis 2014) through the CIPRES Science Gateway v.3.3 interface (Miller et al. 2010). In each Bayesian analysis we specified the following parameters: mcmcp ngen=60000000, burninfrac=0.25, printfreq=6000, and samplefreq=6000; while in each Maximum Likelihood analysis we specified the GTRGAMMA model of nucleotide substitution and 1000 bootstrap iterations. We included Sceloporus grammicus as the sister group to the entire torquatus species group, S. mucronatus as a member of the torquatus species group, as well as Sceloporus bulleri as the sister species of S. torquatus (Martínez-Méndez and Méndez-De la Cruz 2007).

We used Tracer v.1.7.1. (Rambaut et al. 2018) to check the Markov chains (MCMC) convergence implemented in MrBayes, and FigTree v1.4.4 (Rambaut 2018) to visualize the resulting phylogenetic trees.

Genetic distances were calculated using the concatenated matrix of mtDNA data (12S + ND4). Using MEGA X v.10.0.5 (Kumar et al. 2018), we constructed a Neighbor-Joining tree with 1000 bootstrap iterations and the Kimura2-parameter model (Kimura 1980) to subsequently calculate the genetic distances between groups defined by lineages recovered in the phylogenetic analyses under the same parameters. We designed S. bulleri as the external group because it is the sister species of S. torquatus.

We estimated divergence times between lineages using BEAST v2.5.1 (Bouckaert et al. 2019) under a Yule tree model. We inferred models of substitution and rate heterogeneity using bModelTest (Bouckaert and Drummond 2017) for four partitions: 12S, the ND4 coding region, the noncoding tNRA region of ND4, and RAG1. We estimated two separate uncorrelated relaxed clock models for the combined mitochondrial loci and the nuclear RAG1 loci, respectively. A secondary calibration was used to calibrate the node corresponding to the most recent common ancestor between the torquatus species group and S. grammicus. A uniform prior between 12.9 and 18 mya was used for this node, as this range encompasses the estimated divergence date for these taxa in two previous studies on the group (Wiens et al. 2013; Leaché et al. 2016), and has been used in recent divergence estimations for the torquatus species group (Lambert et al. 2019). Three independent runs of 40000000 MCMC generations were run, sampling every 4000 generations. We assessed convergence in Tracer v.1.7.1 (Rambaut et al. 2018), where we compared replicate runs for similar parameter values and then combined them using LogCombiner after discarding the first 10% of trees of each run as burn-in. We used TreeAnnotator to create a maximum clade credibility tree using the median ancestor height and visualized the resulting tree in FigTree v1.4.4 (Rambaut 2018).

We performed a series of statistical analyses to evaluate the multivariate niche overlap between lineages in the environmental spaces. We used the “PCA-env” approach (Broennimann et al. 2012) implemented in the ecospat R packages (Di Cola et al. 2017). This approach calculates niche overlap using the Schoener’s D metric from the first two principal component analysis (PCA) including climate information from the respective lineage occurrence distributions and their background from the calibration area (see below). A smoothed occurrence density was estimated for each lineage using a kernel density function, and this was used to calculate niche overlap. We implemented randomization tests to assess niche similarity for each lineage pair (Di Cola et al. 2017). Here we test whether lineage pairs are more similar than expected based on their background environments (i.e., species are occupying niches that are more similar given the environmental availability in the region). As we are interested in testing a scenario of ecological niche conservatism, testing whether lineages in the S. torquatus complex were more similar than expected by the background conditions is the most appropriate null hypothesis here. We used 100 random replications for these tests. We used an ensemble approach given the high uncertainty in model algorithm selection on transferability under past climate change scenarios. We selected a set of bioclimatic variables for model fit based on collinearity, which was calculated using the Variance Inflation Factor (VIF; Marquaridt 1970). The VIF was calculated for the 19 bioclimatic variables from WorldClim using the vifcor R function from the usdm package (Naimi et al. 2014). Afterward, we selected the following variables for model fit: bio4, bio9, bio15, bio18 and bio19. We evaluated our models by creating pseudo-absences and with data-splitting methods. First, we randomly partitioned the presence data into two sets for calibration (70%) and validation (30%). For each dataset (calibration and validation), we generated a set of pseudo-absences using the ecospat.rand.pseudoabsences function from ecospat R package (Di Cola et al. 2017). The number of pseudo-absences for calibration was 10 times the number of training presences and for validation was 100 times the number of testing presences (i.e., 800 pseudo-absences). Pseudo-absences were created randomly across the entire calibration area or accessible area (M area; Soberón and Peterson 2005) with a minimum distance of at least 5km with respect to presence records. This area represents the hypothetical historical suitable area (HSA) where lineages recovered in our phylogenetic analysis evolved through time. We adopted this validation approach to maximize the number of pseudo-absences in both cross-validation splits and external validation. We used eight model algorithms available in the sdm R package (Naimi and Aráujo 2016), including MaxEnt (Maximum Entropy), MARS (Multivariate Adaptive Regression Splines), GBM (Gradient Boosting Machine), RF (Random Forest), CART (Classification and Regression Trees), SVM (Support Vector Machines), GLM (Generalized Linear Model) and GAM (Generalized Additive Model). Models were trained using 5-folds of cross-validation and 10 bootstrapping replications for a total of 50 replications per algorithm. For each individual model, we evaluated geographical predictive accuracy using the true skill statistic (TSS) and omission rate (Allouche et al. 2006; Fielding and Bell 1997). Finally, we generated a consensus ensemble model weighting for those models maximizing TSS values. This model identifies areas where those models with the highest predictive capacity tend to agree with the environmental conditions for successful population establishment (i.e., habitat suitability distribution).

Then, ensembles were transferred to past climate change scenarios from the paleoclimatic database PaleoClim (Brown et al. 2018) to generate past suitable conditions. This database contains bioclimatic information for 11 time horizons since the last Meghalayan until the mid-Miocene. The time periods are as follows (in parentheses the estimated time period): Meghalayan (4.2–0.3 kya), Northgrippian (8.3–4.2 kya), Greenlandian (11.7–8.3 kya), Younger Dryas Stadial (12.9–11.7 kya), Bølling-Allerød (14.7–12.9 kya), Heinrich Stadial (17.0–14.7 kya), Last Glacial Maximum (LGM; ~21 kya), Inter-Glacial (LIG; ~121 kya), the Marine Isotope Stage 19 in the Pleistocene (MIS19; ~787 kya), mid-Pliocene Warm (~3.2 mya) and the Marine Isotope Stable in the Late Pliocene (M2; ~3.3 mya). These periods include several abrupt global climate change events (Thornalley et al. 2010; 2011; 2013; Brown et al. 2018). We stacked individual models and then estimated the median of suitability values across the region to identify areas where the optimal niche conditions coincided for the majority of lineages as the historical stable areas (HSA).

We tested whether those lineages recovered by molecular phylogenetic analyses exhibit morphological differences through PCA and nMDS methods using morphometric and scutellation characters.

We followed Olson (1990) and Smith (1939) for morphometric and scutellation terminology. All measurements and counts were made by the same person (GCG) using a Mitutoyo 500-196-30 digital caliper (with an accuracy of ± 0.1mm), a 3× magnifier, and a Zeiss 5× stereomicroscope.

We measured 576 individuals exceeding 70mm SVL snout-vent length as S. torquatus reaches sexual maturity at this body size (Guillette and Méndez-De la Cruz 1993; Feria Ortiz et al. 2001). We built a data matrix with 10 morphometric traits (See supplementary file 3: Morphometric measurements). Additionally, a data matrix with 18 scutellation characters (See supplementary file 4: Scutellation counts) was built from 638 adult and juvenile specimens, as these traits are not body size dependent.

We removed the effect of body size on morphometric variables following Velasco and Herrel (2007) where each variable was log10-transformed and regressed against snout-vent length (log10-transformed). The residuals of all variables and the snout-vent length (log10-transformed) were used in a PCA. Then, we performed a Multivariate Analysis of Variance (MANOVA) with the scores obtained from the principal components (PC) to test for significant differences (p <0.05) between means of the variances of the lineages compared.

Alternatively, with the scutellation data matrix we implemented a non-Metric Multidimensional Scaling (nMDS) analysis with the Manhattan coefficient to calculate total differences of the measured variables between individuals of each recovered lineage.

We carried out these statistical analyzes with the tools provided in PAST v.4.01 (Hammer et al. 2001).

We discovered that more than one species is represented in the type series of S. torquatus (Fig. 2). Specifically, the specimen ZMB 628 has divided supraocular scales, 32 dorsal scales, 43 ventral scales, and blue coloration on the belly, throat, and both sides of the head; furthermore, dorsal scales are bordered with black, and light borders of the dark nuchal collar are complete. These characters led us to re-determine this specimen as Sceloporus aureolus Smith, 1942.

Additionally, we re-determined the specimen ZMB 630, a syntype of S. torquatus, as S. t. melanogaster by having undivided supraocular scales, 30 dorsal scales, 41 ventral scales, diffuse dark nuchal collar interrupted by dorsolateral light bands or marks, as well as a series of dark irregular spots that fade over the base of the tail.

Finally, we found that specimen ENCB 5756, a paratype of S. t. mikeprestoni, actually pertains to Sceloporus minor Cope, 1885. This specimen has divided supraocular scales, 36 dorsal scales, 40 scales around the body, and 44 ventral scales.

We obtained 170 sequences from the 12S (321–351 bp), ND4 + adjacent tRNA (553–719 bp), and RAG1 (909 bp) regions. The mitochondrial data matrix contains 60 samples, 1070 bp, 770 conserved sites, 300 variable sites, and 196 parsimony informative sites, while the combined data matrix contains 50 individuals, 1979 bp, with 1639 conserved sites, 350 variable sites, and 205 parsimony informative sites.

The optimal partitioning schemes of the mitochondrial and combined data sets, as well as the best substitution model for each partition, are shown in Table 3.

Mitochondrial gene trees resulting from the Bayesian and ML analyses maintain a similar topology (Fig. 3). We can identify eight different lineages comprising the S. torquatus complex: S. t. torquatus (Posterior Probability, PP=1; Bootstrap, BS=89), S. t. melanogaster (PP=0.99; BS=60), S. t. binocularis (PP=1; BS=100), S. t. mikeprestoni (PP=1; BS=100), S. t. madrensis north (PP=1; BS=100), S. t. madrensis south (PP=1, BS=97), Sceloporus sp. (PP=1, BS=100) and Sceloporus sp. Zacatecas (PP=1, BS=97). The S. torquatus complex was found to be monophyletic with respect to the included outgroup taxa, although with low support (PP=0.77; BS=46).

In both mitochondrial trees, S. t. torquatus, S. t. binocularis, S. t. mikeprestoni, S. t. madrensis north, and S. t. madrensis south forms a clade sister to the clade including Sceloporus sp. and Sceloporus sp. Zacatecas.

Combined mitochondrial and nuclear data phylogenies (Fig. 4) recovered the same eight lineages: S. t. torquatus (PP=1; BS=97), S. t. melanogaster (PP=0.97; BS=32), S. t. binocularis (PP=1; BS=100), S. t. mikeprestoni (PP=1; BS=100), S. t. madrensis north (PP=1; BS=100), S. t. madrensis south (PP=1, BS=100), Sceloporus sp. (PP=1, BS=100) and Sceloporus sp. Zacatecas (PP=1, BS=97). In the ML phylogeny, Sceloporus sp. and Sceloporus sp. Zacatecas form the sister clade to S. t. melanogaster.

We consistently recovered S. bulleri and S. mucronatus as the sister species of the S. torquatus complex, while S. grammicus is sister to all of them.

The genetic distance between S. bulleri and any member of the S. torquatus complex ranges from 0.069–0.085. The genetic distance between S. t. torquatus and S. t. melanogaster is 0.054, between S. t. binocularis and S. t. mikeprestoni is 0.025, between S. t. madrensis north and S. t. madrensis south is 0.045, and that between Sceloporus sp. and Sceloporus sp. Zacatecas is 0.032 (Table 4).

The BEAST time-tree recovered a similar topology and support values to the RAxML and MrBayes trees (Fig. 5). The crown age for the S. torquatus complex is ~5.51 mya (3.61–7.77, 95% HPD). The split between S. t. melanogaster and the two lineages from the Sierra Madre Occidental, Sceloporus sp. and Sceloporus sp. Zacatecas, dates to ~4.33 mya (2.63–6.32, 95% HPD). The divergence between S. t. torquatus and the clade including the four lineages from the Sierra Madre Oriental is ~4.12 mya (2.64–6.03, 95% HPD). The four lineages from the Sierra Madre Oriental, S. t. madrensis north, S. t. madrensis sur, S. t. binocularis, and S. t. mikeprestoni, are recovered as monophyletic with good support and a crown age of ~3.31 mya (2.06–4.76, 95% HPD); relative splitting of the two S. t. madrensis lineages is uncertain, given the low internal posterior probability value within this subclade. The timing of these divergences from the binocularis subclade (S. t. binocularis + S. t. mikeprestoni) is recovered between 1.89–4.4 mya, and the most recent common ancestor between S. t. binocularis and S. t. mikeprestoni lineages is recovered at ~1.81 mya (0.97–2.81, 95% HPD).

In general, there are no similarities in ecological niches of each lineage within the S. torquatus complex (Table 5). Comparison between S. t. melanogaster and S. t. binocularis shows the highest niche similarity (Schoener’s D=0.44), although their respective p-values in the randomization test are discrepant (p=0.01, p=0.09), and therefore this similarity must be taken with reservations.

Potential distribution models (Fig. 6) illuminate some interesting patterns. For example, S. t. torquatus has a greater affinity with existing climatic conditions of central and southern Mexico. Given the suitability values observed in each model, there seems to be reciprocity between the potential areas of S. t. binocularis, S. t. mikeprestoni, and S. t. madrensis north, which together inhabit northeastern Mexico, with respect to the potential area of S. t. madrensis south which is distributed in central eastern Mexico.

According to models projected into the past, the HSA have been very dynamic as they have expanded and contracted consecutively since the late Pliocene, but have remained associated with the main mountainous regions of central and northern Mexico. Between the mid-Pliocene Warm period (~3.2 mya) and MIS19 (~787 kya) another HSA appears in Northeast Mexico. Through the different temporal scenarios, except for the mid-Pliocene Warm, an extensive HSA has been maintained in central Mexico (Fig. 7).

A summary of the descriptive statistics for each taxa is shown in supplementary file 5: Geographic distribution, morphometrics and scutellation of the S. torquatus complex.

We performed a PCA and nMDS analyses to contrast the morphology of the S. torquatus complex members. There is not clear segregation of the analyzed datasets (Figs 8–9).

In the PCA (Fig. 8), the two first PCs explain 70.8% of the total variance, and the MANOVA performed with the scores of all ten PCs yielded Wilks λ=0.4175, F=7.541, p <0.05. The scores, eigenvalues, and percentage of the explained variance are shown in the Supplementary file 6: PCA statistics.

The nMDS analysis (Fig. 9) yielded the following values: Stress value = 1.572; Coefficients of determination (R²): Axis 1 = 0.101, Axis 2 = 0.07448.

According to the International Code of Zoological Nomenclature (ICZN; The International Trust for Zoological Nomenclature 1999) the fixation of a type specimen serves as an objective reference for the application of the taxonomic name it carries (Art. 61.1), and such objectivity is hierarchically continuous from the species level to the family level (Art. 61.1.2). Now, if in the original description of a nominal taxon a specimen or specimens bearing the name was not designated, it is possible that such a designation was made later by the figure of the first reviewer (Arts. 24.2.1). In this context, Taylor (1969) served as the first reviewer designating four syntypes for S. torquatus. The results we present here, show that the specimens ZMB 628 and ZMB 630 belong to distinct taxonomic species other than S. torquatus, thus causing instability in the application of the species name, and therefore warranting a lectotype to be designated from the syntypes (ICZN Arts. 70.3, 74.1). For this purpose, we designate as the lectotype for the name Sceloporus torquatus Wiegmann, 1828 the specimen ZMB 629, and as the paralectotype the specimen ZMB 631. We base the designation of the lectotype on its similarity to the specimen illustrated in Wiegmann (1834; tab. VII, fig. 1), according to ICZN Art.72.4.1.1.

In the other case, misidentification of the S. t. mikeprestoni paratype ENCB 5756 in the original description (Smith and Álvarez 1974) does not exclude it from the type series of this nominal taxon (ICZN Art.72.4.2).

This study includes genetic data from S. t. mikeprestoni and S. t. madrensis for the first time ever, as well as the most extensive sampling throughout the distribution of the S. torquatus complex, to accomplish the most complete molecular phylogeny of this emblematic group of phrynosomatid lizards to date. The S. torquatus complex is a monophyletic group composed of eight independent lineages, five of which represent recognized subspecies, while the remaining three represent unnamed taxa that are awaiting descriptions (Flores-Villela et al. in prep.).

There is evidence to recognize the southern populations of S. t. madrensis as an independent lineage, previously confused with S. t melanogaster and S. t. madrensis (Smith 1939; Olson 1991). Therefore, there is no zone of sympatry between the southern populations of S. t. madrensis and S. t. torquatus, as Olson (1991) argued.

In addition, we confirm the existence of another cryptic species from western Mexico suggested by Martínez-Méndez and Méndez-De la Cruz (2007) and Martínez-Méndez et al. (2019). In our research, the lineages Sceloporus sp. and Sceloporus sp. Zacatecas are more closely related to S. t. melanogaster given the smaller genetic distances between them and their geographical proximity (Figs 1, 5; Table 4).

Additional tissue samples from the northernmost populations of S. t. melanogaster could help elucidate phylogenetic relationships within the S. torquatus complex as a sister taxon of Sceloporus spp. from Nayarit, Jalisco and Zacatecas, and could also solve phylogenetic relationships within the lineage S. t. melanogaster. Future samplings along the contact zone of S. t. torquatus and S. t. melanogaster in central Mexico would be useful to determine the extent of gene flow, and to investigate mechanisms of reproductive isolation, especially since behavior and coloration are known to be related to conspecific recognition and reproductive success in Sceloporus (Hunsaker 1962; Jiménez-Arcos et al. 2017).

At the end of the Neogene, tectonic and volcanic activity gave rise to the main mountain systems of Mexico, promoting vicariant events in numerous taxa (Morafka 1977; Bryson et al. 2012). Our estimation of divergence times shows that the current phylogeographic structure of the S. torquatus complex coincides with this period. Pleistocene climate changes may have led to the diversification of numerous taxa (Bryson et al. 2011, 2012; Leaché et al. 2013; Díaz-Cárdenas et al. 2019), to which we include the most recently diverged lineages within the S. torquatus complex, S. t. binocularis and S. t. mikeprestoni. Geographically, the nearest localities for S. t. binocularis and S. t. mikeprestoni are separated by ~30km airline and ~1600m in elevation.

The modeled HSA (Fig. 7) indicates that climatic conditions have been favorable for the S. torquatus complex repeatedly in east-central Mexico, as far east as Veracruz, very close to the Gulf of Mexico slope. We did not find populations of S. torquatus beyond eastern Tlaxcala during field work, and the historical records from Veracruz that we examined were redetermined as S. mucronatus and S. formosus. According to the HSA, future sampling in western Mexico, along the Sierra Madre Occidental could reveal the discovery of new species of Sceloporus related to the S. torquatus complex, as was the case of Sceloporus spp. from Nayarit, Jalisco and Zacatecas.

Although the niche similarity test that we performed is not conclusive (Table 5), we note the fact that four main biogeographic provinces constitute the current geographic distribution of the S. torquatus complex, implying a great heterogeneity of suitable habitats and topography. The distribution of some lineages within the S. torquatus complex are isolated in mountainous areas where the climatic change has been accelerated (Sinervo et al. 2010), thus future studies could be directed to reassess the extinction risk of these lineages.

It has been suggested that morphological convergence may be related to environmental similarity in other species of the torquatus group (Martínez-Méndez et al. 2012). In the case of the S. torquatus complex, we observed morphological conservatism (Figs 8–9) despite their wide geographic distributions (Fig. 1) and heterogeneity of environments inhabited (Table 5). This is clearly observed in S. t. torquatus which is distributed in central Mexico, where it inhabits mainly pine forests, oak forests and scrub, at 1300–3533m, and S. t. melanogaster which is distributed throughout central and much of northern Mexico, living mainly in different types of scrub, pine forests, oak forests and grasslands, at 1100–2745m. Both are found frequently on stone walls or fences, in agricultural land, and urban areas. The lack of considerable differences in morphometry and scutellation among the populations of the compared lineages may be a consequence of the relatively recent diversification of the S. torquatus complex. In live and preserved specimens, coloration characteristics are generally useful for distinguishing between members of the S. torquatus complex, except for specimens from the wide contact zone between S. t. torquatus and S. t. melanogaster along the Faja Volcánica Transmexicana.

As we expected due to its inherent properties, the mtDNA showed higher genetic differentiation than the nDNA and largely drove the phylogeographic patterns discussed above. While we acknowledge the limitations of using solely or mainly mtDNA for species delimitation (Leaché and Mulcahy 2007; Leaché 2010), this practice is commonly used to discern recently diverged lineages such as those comprising species complexes, even within Sceloporus (Wiens and Penkrot 2002; Martínez-Méndez et al. 2012; Díaz-Cárdenas et al. 2017), that would otherwise be difficult to resolve with nDNA alone. The inclusion of genome-wide nuclear markers, such as those generated from next-generation sequencing technology, would add further phylogenetic and taxonomic resolution to this group.