Taxonomic hypotheses and the biogeography of speciation in the Tiger Whiptail complex (Aspidoscelis tigris: Squamata, Teiidae)

(225)Biodiversity in southwestern North America has a complex biogeographic history involving tectonism interspersed with climatic fluctuations. This yields a contemporary pattern replete with historic idiosyncrasies often difficult to interpret when viewed from through the lens of modern ecology. The Aspidoscelis tigris (Tiger Whiptail) complex (Squamata: Teiidae) is one such group in which potential taxonomic boundaries have been confounded by a series of complex biogeographic processes that have defined the evolution of the clade. To clarify this situation, we first generated multiple taxonomic hypotheses, which were subsequently tested using mitochondrial DNA sequences (ATPase 8 and 6) evaluated across 239 individuals representing five continental members of this complex. To do so, we evaluated the manner by which our models parsed phylogenetic and biogeographic patterns. We found considerable variation among species ‘hypotheses’, which we suggest in part reflects inflated levels of inter-population genetic divergence caused by historical demographic expansion and contraction cycles. Inter-specific boundaries with A. marmoratus juxtaposed topographically with the Cochise Filter Barrier that separates Sonoran and Chihuahuan deserts (interpreted herein as case of ‘soft’ allopatry). Patterns of genetic divergence were consistent across the Cochise Filter Barrier, regardless of sample proximity. Surprisingly, this also held true for intraspecific comparisons that spanned the Colorado River. These in turn suggest geomorphic processes as a driver of speciation in the tigris complex, with intraspecific units governed locally by demographic processes. HIGHLIGHTS Phylogeographies of vertebrates within the southwestern deserts of North America have been shaped by climatic fluctuations imbedded within broad geomorphic processes. The resulting synergism drives evolutionary processes, such as an expansion of within-species genetic divergence over time. Taxonomic inflation often results (i.e., an increase in recognized taxa due to arbitrary delineations), such as when morphological divergences fail to juxtapose with biogeographic hypotheses. However, isolated groups can be discriminated within-species by mapping genetic variability onto geographic distances. This approach can often diagnose ‘hard’ barriers to dispersal, or alternatively, strong selection acting against hybridization. On the other hand, elevated genetic divergences among groups less-isolated would underscore isolation-by-distance (i.e., an increase in genetic differentiation concomitant with geographic distance). The biogeographic patterns we identified in Tiger Whiptail are largely synonymous with those found in other regional species, particularly given the geomorphic separation of Mohave and Sonoran deserts by the Colorado River, and Sonoran/ Chihuahuan deserts by the Cochise Filter Barrier. Our results for the Tiger Whiptail complex broaden and extend the context within which polytypic species are conserved and managed, particularly those that reflect an incongruence among molecular and morphological standards.

(225)Biodiversity in southwestern North America has a complex biogeographic history involving tectonism interspersed with climatic fluctuations. This yields a contemporary pattern replete with 2 historic idiosyncrasies often difficult to interpret when viewed from through the lens of modern ecology. The Aspidoscelis tigris (Tiger Whiptail) complex (Squamata: Teiidae) is one such 4 group in which potential taxonomic boundaries have been confounded by a series of complex biogeographic processes that have defined the evolution of the clade. To clarify this situation, we 6 first generated multiple taxonomic hypotheses, which were subsequently tested using mitochondrial DNA sequences (ATPase 8 and 6) evaluated across 239 individuals representing 8 five continental members of this complex. To do so, we evaluated the manner by which our models parsed phylogenetic and biogeographic patterns. We found considerable variation among 10 species 'hypotheses', which we suggest in part reflects inflated levels of inter-population genetic divergence caused by historical demographic expansion and contraction cycles. Inter-specific 12 boundaries with A. marmoratus juxtaposed topographically with the Cochise Filter Barrier that separates Sonoran and Chihuahuan deserts (interpreted herein as case of 'soft' allopatry). 14 Patterns of genetic divergence were consistent across the Cochise Filter Barrier, regardless of sample proximity. Surprisingly, this also held true for intraspecific comparisons that spanned the 16 Colorado River. These in turn suggest geomorphic processes as a driver of speciation in the tigris complex, with intraspecific units governed locally by demographic processes. 18 HIGHLIGHTS(185) 20 1. Phylogeographies of vertebrates within the southwestern deserts of North America have been shaped by climatic fluctuations imbedded within broad geomorphic processes. 22 2. The resulting synergism drives evolutionary processes, such as an expansion of withinspecies genetic divergence over time. Taxonomic inflation often results (i.e., an increase 24 in recognized taxa due to arbitrary delineations), such as when morphological divergences fail to juxtapose with biogeographic hypotheses. 26 3. However, isolated groups can be discriminated within-species by mapping genetic variability onto geographic distances. This approach can often diagnose 'hard' barriers to 28 dispersal, or alternatively, strong selection acting against hybridization. On the other hand, elevated genetic divergences among groups less-isolated would underscore 30 isolation-by-distance (i.e., an increase in genetic differentiation concomitant with geographic distance). 32 4. The biogeographic patterns we identified in Tiger Whiptail are largely synonymous with those found in other regional species, particularly given the geomorphic separation of 34 Mohave and Sonoran deserts by the Colorado River, and Sonoran/ Chihuahuan deserts by the Cochise Filter Barrier. 36 5. Our results for the Tiger Whiptail complex broaden and extend the context within which polytypic species are conserved and managed, particularly those that reflect an 38 incongruence among molecular and morphological standards.  (Barley et al. 2019). However, intraspecific morphological variation coupled with high rates of hybridization has made species-delineation difficult, which in turn has 48 muddied their taxonomy. The occurrence of paraphyly in the group was partially resolved by subsequently partitioning species between a resurrected Aspidoscelis and the nominate 50 Cnemidophorus (Reeder et al., 2002). While a positive step, it also necessitated reallocation of several groups, such as: Aspidoscelis cozumelus (parthenogenetic species), A. 52 deppi (gonochoristic species), A. sexlineatus (parthenogenetic and gonochoristic species), A. tesselatus (parthenogenetic species), and A. tigris (gonochoristic species). As a result, many 54 systematic issues inherent within Cnemidophorus were simply transferred to Aspidoscelis.
The A. tigris species-group, for example, represents a perplexing systematic challenge as 56 it comprises a widely distributed complex of gonochoristic populations located within desert regions of southwestern North America, northern México (to include Baja California), as well as 58 numerous islands east and west of the peninsula. Although the group per se is characterized by a readily assessed set of morphological characters (Walker and Maslin, 1981) and a distinctive 60 karyotype (Lowe et al., 1970), its diversity (i.e., ranging from a single species to more than a dozen) has been a contentious issue in North American Herpetology. Much of its diversity is 62 represented by distinct populations located on relatively small, vulnerable islands in the Gulf of California. Several of these were designated as distinct species (Walker et al., 1966a, b;Walker 64 and Taylor, 1968;Walker and Maslin, 1981), based snout-vent length, color pattern, and/or morphology (e.g., A. bacatus of Isla San Pedro Nolasco; A. catalinensis of Isla Santa Catalina; A. 66 martyris of Isla San Pedro Martír; and A. celeripes of Isla San Jose). However, some authors (e.g. Lowe et al., 1970;Wright, 1993) have considered them instead as subspecific variants. 68 Continental populations have also proven taxonomically contentious. For example, a longstanding disagreement (Hendricks and Dixon, 1986;Dessauer et al., 2000) centers on whether A. 70 tigris and A. marmoratus in small areas in Hidalgo County NM represent intergrading subspecies or hybridizing species. 72 Taxonomic difficulties are thus quite apparent within the A. tigris complex, and primarily stem from complex patterns of intergradation and hybridization confounded by non-74 discriminating morphological variation. To provide clarification, we evaluated a complex of four continental subspecies (A. t. mundus,A. t. punctilinealis,A. t. septentrionalis,and A. t. tigris;76 Fig. 1), plus a recently elevated putative sister species, A. marmoratus (Hendricks and Dickson, 1986;Reeder et al., 2002;Tucker et al., 2016). Our goal was to examine how mitochondrial 78 (mt)DNA variation in A. tigris subspecies compares with morphological discordance [e.g., subspecies seemingly supported by color pattern yet weakly differentiated by scalation and 80 morphometrics (Taylor et al., 1994;Walker et al., 2015;Fig. 2)]. Furthermore, morphological intergrades and overlapping genetic polymorphisms suggest admixture among several groups: A. 82 t. tigris and A. t. septentrionalis (Taylor 1988), A. t. tigris and A. t. punctilinealis (Taylor 1990), and A. t. punctilinealis and A. t. septentrionalis (Dessauer et al. 2000). 84 The contact zone separating A. t. punctilinealis and A. marmoratus (i.e., A. t. marmoratus of the following authors) coincides with a genetic cline at multiple loci Cole 1991, 86 Dessauer et al. 2000), suggesting secondary contact following a period of isolation (Barton and Hewitt, 1985). A steep cline involving six color-pattern characteristics is also apparent in the 88 hybrid zone between A. t. punctilinealis and A. marmoratus (Dessauer et al., 2000), arguing for the utility of these characters for delineation of A. tigris lineages. 90 In the present study, we focused on defining phylogeographic patterns and interpreting phylogenetic delineations primarily within the range of A. t. punctilinealis in Arizona. To the 92 east, it forms steep morphological/ genetic clines at a hybrid zone with A. marmoratus (Dessauer et al., 2000), whereas to the west it exhibits a relatively gradual phenotypic gradient with A. t. 94 tigris (Taylor 1990). Interestingly, both A. t. septentrionalis and A. t. punctilinealis also form exceptionally gradual morphoclines with A. t. tigris that lack apparent ecological underpinnings 96 (Taylor et al. 1994, Taylor 1990). This stands in contrast to that found between A. t. punctilinealis and A. marmoratus (<2km;Dessauer et al., 2000;Zweifel 1962). Given this 98 discontinuity, we also examined the association between morphological gradation in the complex and its genetic discontinuity. To do so, we extended the phylogenetic resolution of lineages 100 within the delineated A. tigris complex, then contrasted our phylogenetic patterns with morphological species identifications. We generated taxonomic hypotheses by employing 102 'species delimitation' algorithms, then employed subsequent downstream analyses to compare the manner by which each model partitioned genetic variation. 104 The Colorado Plateau as a biogeographic entity 106 The Colorado Plateau is a tectonically intact landform (~2000m elevation) with scant internal deformation, thus providing a relatively linear configuration of layered sedimentary rock 108 for the Plateau interior that contrasts with an eroding margin since Late Miocene (Levander et al., 2011). It also stands in sharp contrast with the more fractured geomorphic features of 110 neighboring regions: The Rocky Mountains (immediately north and east), Rio Grande Uplift (southeast), and the Basin and Range Province (immediately west) ( Fig. 1 of Flowers, 2010). The 112 uplift of the Plateau and the integration of the Colorado River were defining geomorphic events (Minckley et al., 1986;Spencer et al., 2008). 114 Prior to integration, the Colorado River emptied into an endorheic basin, most likely following the early path of a major modern tributary (Little Colorado River; House et al., 2008). 116 In late Miocene-early Pliocene, headwater erosion and basin 'spillover' forged a connection between the river and the Gulf of California via the Salton Trough (Dorsey et al., 2018;118 Nicholson et al., 2019;Sarna-Wojcicki et al., 2011). The Gila River system, which had previously drained directly into the Gulf (Eberly & Stanley, 1978) was similarly incorporated. 120 Other drainages that terminated into closed basins during the Miocene (i.e. Salt and Verde rivers) were also integrated with the Gila River system by late Pliocene (Blakely & Ranney, 2008). The 122 Pliocene also brought decreased seasonality as well as elevated temperatures and humidity, a climatic regime subsequently reversed by Quaternary glaciation (Thompson, 1991). 124 Geomorphic events were not only pivotal with regard to biological diversification on the Plateau (Bangs et al., 2020a(Bangs et al., , 2020bDouglas et al., 1999;Douglas et al., 2016) but also 126 provoked endemism as well, in synergy with widespread physiographic and climatic oscillations.

Phylogenetic inference 158
To develop initial hypotheses, we first inferred phylogenetic relationships using both Bayesian and distance-based methods. We tested models of nucleotide substitution (N=88-160 JMODELTEST; Darriba et al., 2012), using Akaike and Bayesian information criteria (AIC and BIC, respectively). We generated an ultrametric tree using the HKY-model (selected by both 162 AIC and BIC) with gamma distributed rates and invariant sites (BEAST2; Bouckaert et al., 2014). We applied a Yule model with default parameters and a relaxed log-normal clock process. 164 MCMC chains were sampled every 500 th of 10 million MCMC iterations, following a burn-in period of 2 million. Chain convergence and mixing were visually assessed, and MCMC samples 166 thinned so as to maintain effective sample sizes >200 while reducing autocorrelation of samples (TRACER; Rambaut and Drummond, 2007). The posterior sampled trees were summarized as a 168 maximum clade credibility tree with branch lengths representing median heights. For contrast, we also generated a neighbor-joining (NJ) tree using model-corrected distances (R-package APE; 170 Paradis et al., 2004) in R (R Core Team, 2020), with nodal support assessed via 100 bootstrap pseudo-replicates, as well as a maximum-likelihood tree in PHYML (Guindon et al. 2010). Basic 172 population summary statistics for major phylogenetic clades (as diagnosed via nodal support and statistical tests of monophyly) were also computed in PEGAS (Paradis, 2010). 174

Delimiting groups within A. tigris 176
We first evaluated the currently accepted taxonomic designations for subspecies within A.
tigris (e.g., Crother et al., 2017; hereby abbreviated as 'TAX'), then employed an array of 178 analytical approaches to generate six additional delimitation hypotheses. Distance-based clustering methods, used in the context of species delimitation, assume that biological units (i.e., 180 species) form non-overlapping genetic clusters. We first employed a recursive method for 'Automated Barcode Gap Discovery' (ABGD; Puillandre et al., 2012) that statistically infers a 182 threshold distance from the data (i.e., the genetic distance at which groups are considered to represent different clusters). We also generated clusters using corrected genetic distances as input 184 (CD-HIT; Fu et al., 2012), with manually specified distance thresholds of 2% and 5% (hereafter referred to as CD2 and CD5). These hard thresholds were chosen because they represent 186 approximate 'breakpoints' at which changes occur in the ABGD analysis regarding the number of sampled groups. 188 We defined phylogenetic species using several approaches: The first, a 'naïve' method, defined groups as the shallowest divergences statistically supported and reciprocally 190 monophyletic (Rosenberg, 2007), hereafter abbreviated as 'ROS' ( SPIDER R package; Brown et al., 2012). A second approach derived groups from tree characteristics under an assumed 192 speciation process. The first of these methods ['Generalized Mixed Yule Coalescent' (GMYC; Pons et al., 2006)] represents the species boundary as a distinguishable change in branching 194 patterns within the phylogeny, with independently evolving clusters diagnosed as those at which branching patterns transition from expected under two different models: An assumed 196 interspecific (Yule process; Nee et al., 1994) and an intraspecific (neutral coalescent; Hudson, 1990). In the second ['Poisson Tree Process' (PTP; Zhang et al., 2013)] the probability of 198 speciation is considered in terms of mutational distance between cladogenesis events (as assumed with a Poisson distribution). PTP requires a tree with branch lengths scaled per 200 substitutions, and given this, our input represented a tree topologically consistent with that inferred by BEAST but derived instead with PHYML (Guindon et al. 2010). GMYC was run 202 using the Splits R package (Fujisawa et al., 2013); PTP using the Bayesian implementation provided by the bPTP server (Zhang et al., 2013). 204

Comparingf delimitations 206
We used 2-way and 3-way Analyses of Molecular Variance (AMOVA; Excoffier et al. 1992) with group assignments defined by individual membership under each hypothesis. We 208 then used isolation-by-distance (IBD) as a null model under the assumption that the entire complex represented a single species. Our justification stemmed from the pervasive nature of 210 IBD with regard to dispersal-limited species (Sexton et al. 2014). In this sense, deviations from the expected IBD pattern are hypothesized as diagnostic of some alternative non-spatial process 212 (e.g. reproductive isolation). Partial Mantel tests (Smouse et al. 1986), employing pairwise genetic distances with geographic distances as a covariate (Douglas and Endler, 1982;Meirmans,  All clades were significantly supported as being reciprocally monophyletic [Rosenberg (2007) test], and all reflected similarly high bootstrap support in the neighbor-joining analysis. The 234 maximum-likelihood analysis resulted in a tree topologically identical, and was thus not displayed. 236 All four subspecies formed a monophyletic group sister to a well-defined A. marmoratus (clade H). Aspidoscelis t. septentrionalis (clade G, excluding mundus samples) demonstrates 238 significantly negative values for both Tajima's and Fu and Li's D statistics, with exceptionally low nucleotide diversity and haplotype diversity much less than that recorded for A. t. 240 punctilinealis and A. t. tigris (Table 1). These suggest a population expansion following a bottleneck (Douglas et al., 2010). Aspidoscelis marmoratus showed similar signatures, although 242 with a much greater differential between haplotype and nucleotide diversities, and a reduced number of segregating sites relative to sample size (Table 1) punctilinealis were subdivided into five major clades with strong nodal support and reciprocal monophyly (clades A-F, Fig. 3 Chihuahuan deserts (Castoe et al. 2007).
All clades within A. t. tigris/punctilinealis generally showed higher numbers of 258 segregating sites relative to nucleotide diversity (negative values for Tajima's D; clades A-E Table 1). Clades A-D also showed a high rate of singleton haplotypes. Although non-significant 260 (Table 1) (Taylor et al. 1994 Fig. 4). GMYC and PTP produced large numbers of groups, 292 with 13 in the former and 16 in the latter, with many consisting of but a single haplotype.

Partitioning of genetic variance
The spatial analysis supported models with fewer species. In the bootstrapped datasets, 296 FCT overlapped among K values (Fig. 5) Fig. 6). Significant associations only 302 occurred with the ROS, TAX, ABGD and CD5 models when permutations were limited to those strata defined by delimitation models (i.e., stratified Mantel). The partial Mantel tests with the 304 ABGD model produced the strongest association (r=0.777), with values generally declining as the number of species increased ( Table 2). Correlation of group membership with genetic 306 distance under all models was greater than in the simple Mantel test for all cases, save for the PTP model (r=0.5754), suggesting that all delimitations better encapsulate spatial genetic 308 variation than one in which all examined samples are collapsed.
Overall, those models partitioning A. tigris subspecies were supported by AMOVA 310 analyses (Supplementary Table S1, S2), with high genetic variance among groups rejecting panmixia. Among-group genetic variance (analogous to FCT) varied between 72-90% and was 312 maximized in models that recognizing greater numbers of partitions (i.e., CD2, GMYC, and PTP), with 89-90% of variance explained. These values generally declined as models became 314 more conservative.
Interpolating the AIS results with a diffusion model (Fig. 7)  Cochise Filter Barrier (CFB), a biogeographic entity that has similarly been implicated as influential in defining phylogeographic relationships among other taxa (e.g., Castoe et al., 2007;344 O'Connell et al., 2017;Morafka, 1977). Pyron et al. (2015) invoked this barrier as a potential driver of so-called 'soft allopatry' (Hickerson and Meyer 2008), in which niche differentiation 346 drives speciation in lieu of a hard barrier to dispersal. Such barriers can often lead to permeability at contact zones, given that divergence is maintained by adaptation rather than 348 physical separation. As might be expected, a well-characterized hybrid zone separates A. t. punctilinealis and A. marmoratus, replete with steep clines in color polymorphism (Dessauer et 350 al., 2000). We also identified a distinct east-west divide in mitochondrial haplotypes that is concordant with an hypothesis of soft allopatry. 352 The potential for hybridization 354 We found phylogenetic subdivision less clear among subspecies of A. tigris and surmise this likely reflects a cumulative effect of contemporary physiography coupled with historic 356 processes. Aspidoscelis t. tigris and A. t. punctilinealis were not distinct. Both transverse the Colorado River, and neither is strictly defined by the Mohave-Sonoran transition. This is perhaps 358 not surprising, given that morphological intergradation between A. t. tigris and A. t. punctilinealis does not follow a strict ecological gradient, but instead manifests itself as an 360 exceptionally wide and smooth morphocline (Taylor et al. 1994).
Aspidoscelis t. tigris from the most northerly part of the range had mitochondrial 362 sequences aligning instead with A. t. septentrionalis, potentially suggesting exchange between the two (Taylor, 1988). Here, we emphasize that increased genetic sampling is necessary so as to 364 understand the nature of this putative exchange. In contrast, A. t. punctilinealis and A. t. septentrionalis were demarcated according to recognized physiographic provinces, with A. t. 366 septentrionalis occupying the Colorado Plateau and A. t. punctilinealis the Sonoran Desert (Fig.   4, S1). This was expected, particularly given our understanding of their respective ranges (e.g. 368 Fig. 1) and associations with the disparate thorn-scrub communities therein (e.g., Zweifel, 1962;Pianka, 1970). We suggest this ecological gradient represents a 'soft boundary,' similar to those 370 that mediate the transition between A. t. punctilinealis and A. marmoratus. Unique communities characterize both the Colorado Plateau and Sonoran Desert (e.g., Douglas et al., 2002Douglas et al., , 2010, and 372 this suggests a common, underlying process driving diversification therein (Riddle and Hafner, 2006;Pianka, 1967 A. t. punctilinealis to the southeast. This is likely a result of mitochondrial exchange among these subspecies, again highlighting the need for a genomic evaluation. 384 Taylor et al. (1994) demonstrated the existence of a wide, smooth morphocline >180 km in width that traversed the Colorado River and separated A. t. tigris from A. t. punctilinealis 386 along the southern extent of their contact (Fig. 1A). This stands in contrast to a 16 km-wide morphocline found at the northern extent of their contact, and that discrepancy begs an obvious 388 ecological explanation (Taylor 1990). The northern morphocline coincides with the Arizona-California-Nevada nexus, a hotspot for mitochondrial diversity as well, where sub-clades overlap 390 without clear physical boundaries (Fig. 4, S1). In addition, multiple clades (A and F) traverse the Colorado River, with the highly divergent clade F (found primarily within morphologically 392 defined A. t. tigris) possibly representing a more widespread Mohave clade.
One explanation for the above variation in clinal width evokes an historic perspective 394 wherein the northern contact between A. t. tigris and A. t. punctilinealis is secondary, whereas range overlap to the south was maintained. This would generate morphoclines of varying ages 396 and would offer one potential explanation for the variance in mitochondrial divergence (i.e., elevated to the north, nonexistent to the south), as well as morphocline width. Here, we posit that 398 the duration and extent of contact may be a driving mechanism. A related hypothesis would be that the width of the cline is in response to an underlying environmental gradient, such as the 400 transition between Sonoran and Mohave habitats. However, neither hypothesis can be tested with the data presented herein. 402 An additional possibility is the fluctuating capacity of the Colorado River as a barrier to dispersal. We offer scant evidence for it serving as a vicariant structure, a result corroborated by 404 studies of other reptiles (e.g., desert tortoise: Lamb et al., 1989;iguanids: Lamb et al., 1992).
Although the Colorado River does seemingly constrain dispersal (Pounds and Jackson 1981), 406 gene flow could still be facilitated by the periodic fluctuations in its trajectory and flow. In fact, the shallow divergences found among clades transecting the river provide evidence for how 408 recent these events were.  (Hewitt 2000). This is reflected in an increased number of low-frequency haplotypes, and low nucleotide diversity among 426 contemporary haplotypes (Table 1). We also noted similar patterns in A. t. septentrionalis and A. marmoratus, suggesting bottleneck-expansion events within these taxa as well. 428 The effects of the 'refugium phase' can be seen when examining modern short-range endemics [i.e., those with distributions easily contained within a 100x100 km grid (i.e., <10,000 430 km 2 ]. An example would be the isolated New Mexico Ridge-nosed Rattlesnake (Crotalus willardi obscurus) that reflects signals of relatively recent divergence and population reduction 432 (Holycross and Douglas 2007). Its historic connectivity in the Sky-islands has been truncated by contemporary climate change, whereas its refugium phase, as gauged by demographic and niche 434 modeling, underscored a continued loss of habitat (Davis et al. 2015). We suspect a similar refugium-expansion model would explain the extant genetic differentiation found in the present 436 study.

Comparing species delimitation approaches
We found considerable variance in the number of groups delimited using the various 440 algorithmic methods, especially when contrasted with the current taxonomy (Fig. 3). This could reflect the capacity of the various models to diagnose cryptic lineages, or alternatively, suggest 442 the data violate implicit assumptions (Carstens et al. 2013). Phylogenetic models (i.e., ROS, PTP, GMYC) generally produced the largest numbers of groups, but did not explain the 444 geographic patterns of genetic variation as well as did other models (Table 2; Figure 3). All models, however, captured greater genetic variance among groups (i.e., 72-90% in AMOVA) 446 than within-group (10-28%), indicating the recognition of an authentic biological signal (Table   S1, S2). 448 Species delimitation employing the coalescent, such as GMYC, have previously been criticized for a tendency to over-split when applied to structured populations (Lohse 2009). This, 450 in turn, suggests that a single-threshold delimitation model may be inappropriate for use with hierarchically structured data (Talavera et al. 2013). Strong intraspecific structure in 452 Aspidoscelis, similar to that found in other desert reptiles Douglas 2007, Douglas et al. 2002), may explain the over-fitting of GMYC. A more comprehensive test using additional 454 scenarios is warranted, as this will more appropriately characterize those biases that can emerge when the method is applied to empirical data. 456 We also employed a second phylogenetic method (PTP) that essentially modeled speciation as a clock-like process, with its probability of occurrence represented as a function of 458 time (i.e., branch length). Given this, any aspect that would promote rate heterogeneity could be sufficient to violate the model, such as the extreme demographic fluctuations suspected for 460 Aspidoscelis populations.
Clustering-based methods, on the other hand, may be biased by elevated genetic distances 462 or incomplete sampling, with group delimitation being driven entirely by the arbitrary nature of the (user-selected) threshold. The ABGD model attempts to mitigate this by automatically 464 detecting distance-based partitions, as bounded by priors that account for maximum allowable intraspecific divergence (Puillandre et al. 2012). ABGD and clustering methods appeared to 466 delimit groups that more appropriately represented the genetic/ geographic distances paradigm.
When examining genetic distances defined by these methods within-versus among-group, we 468 found A. marmoratus approximately equal to A. tigris with regard to genetic distances, regardless of whether A. tigris was intact or partitioned (Fig. 6). Among-group genetic distances were equal 470 regardless of how proximate were the samples (i.e., genetic x geographic distance slope approaches zero). This suggests that divergence is maintained not by physical distance but 472 instead by some intrinsic qualities possessed by those groups, such as reproductive isolation or mutually exclusive habitat requirements (Fig. 6). In this case, near-zero slopes in the among-474 group comparisons lend some validity to a species-level designation for A. marmoratus and suggested the potential for reproductive isolation within the tigris complex (i.e., between A. t. 476 tigris/punctilinealis and A. t. septentrionalis; Fig. 6B). However, reduced sample sizes for A. t. mundus and A. t. tigris remain an issue. Within-group comparisons resulted in a positive 478 relationship with geographic distance, as was expected (Fig. 6).

Species concepts and taxonomic conclusions for Aspidoscelis
Species concepts present a philosophical quandary (Coyne and Orr, 2004). For example, 482 many modern interpretations [i.e., biological (BSC; Mayr, 1942); evolutionary (ESC; Frost and Hillis, 1990); phylogenetic (PSC; Cracraft, 1983)] offer a mechanistic or philosophical 484 interpretation for the existence of species, yet do not themselves present an explicit operational criterion that can universally define species. The definitions also have fundamental conflicts with 486 regards to the secondary characteristics each relies upon (e.g., reproductive isolation, morphological diagnosability; de Queiroz, 2007). 488 Given that inferring a direct relationship between genotypic (or 'haplotypic') clusters is a questionable approach, we instead relied upon a 'theory-independent' definition (the 'genotypic 490 cluster' concept: Douglas et al., 2007;Sullivan et al., 2014) wherein clusters are defined as '… distinguishable groups of individuals that have few or no intermediates when in contact' (Mallet, 492 1995). We used this basis to compare among delimitation hypotheses for Aspidoscelis tigris, and how they may (or may not) extend to theoretical species definitions such as the BSC. 494 An expectation of the '… little to no intermediates in contact' criterion is that individuals from separate genotypic clusters should not be any more genetically similar when sampled in 496 proximity than when sampled from more geographically disparate locations. Such a relationship would suggest continuous gene flow limited only by the ability of individuals to disperse (i.e., 498 'isolation-by-distance' sensu Wright, 1940Wright, , 1943. We tested the applicability of this criterion to our species hypotheses using spatial and 500 non-spatial analyses of molecular variance (AMOVA/ SAMOVA), by testing the manner by which each partitioning scheme explained spatial signatures of genetic differentiation (IBD/ 502 Mantel), and with qualitative comparisons of landscape-level genetic differentiation (AIS interpolation). We found the ABGD model (i.e., separation of A. marmoratus from all A. tigris 504 subspecies) did well with partitioning spatial haplotypic variation, with A. marmoratus being highly differentiated (>10% sequence divergence) even from localities sampled <2km from A. t. 506 punctilinealis (Fig. 6A). These patterns are consistent with many commonly applied species concepts (e.g. reproductive isolation; per BSC), and given this, the elevated status of A. 508 marmoratus is sustained.
We also found limited (albeit, existing) exchange between two major clades within the 510 remaining A. tigris: A. t. tigris/punctilinealis versus A. t. septentrionalis/mundus (i.e., TAX model). Here, genetic differentiation was again uncorrelated with spatial distances (Fig. 6B), 512 possibly implicating intrinsic or extrinsic factors that may satisfy species-level criteria. Previous estimates place divergence of these groups to be no greater than 5-10 kya, based upon 514 availability of suitable habitat following the Wisconsin glacial period (Zweifel, 1962). If so, then our observed mean corrected genetic distance of 0.058 equates to a mutation rate of 5.8 516 substitutions per million years, a rate considerably greater than the expected rate of ~0.0045 for ATPase 6, the more rapidly evolving of our two mitochondrial genes (Kumar, 1996;Mueller, 518 2006).
Although we recognize that numerous factors contribute to mtDNA rate heterogeneity, 520 our rudimentary calculation suggests that subspecies within A. tigris may have diverged much earlier than previously thought. Although our results support the continued recognition of A. t. 522 septentrionalis as distinct, we deem them insufficient to propose elevation beyond a subspecific level. Instead, we argue for greater sampling at both geographical and genomic scales to properly 524 adjudicate taxonomies within A. t. tigris. The latter, while enabling tests of hybridization/ introgression within a more statistically robust framework (e.g., Chafin et al., 2019), would also 526 open avenues for novel species delimitation algorithms to be formulated using unsupervised machine learning approaches (Derkarabetian et al., 2019;Martin et al., 2020;Mussmann et al., 528 2020 tigris and A. t. punctilinealis were not reciprocally monophyletic, however assigning causality of 532 this to hybridization/ intergradation, stochastic lineage sorting, or a lack of phylogenetic signal requires additional data (Chafin et al., 2020). Again, we suggest a more in-depth examination 534 using a larger panel of molecular markers so as to confirm or reject continued recognition of A. t. tigris as distinct from A. t. punctilinealis, and likewise to adjudicate recognition of A. t. mundus 536 from A. t. septentrionalis.

Conclusion
We found that algorithmic species delimitation methods vary substantially in the manner 540 by which they partition hierarchical phylogenetic diversity, particularly when contrasted against the respected overlays of historical biogeography, morphology, ecology, and prior genetic 542 evaluations. By this we mean that some algorithms (e.g., ABGD) yielded groupings which support external evidence of divergence at the 'species' level (e.g., according to a BSC-like 544 interpretation), while others (e.g., CD5) yield groupings that could potentially fit various operational species concepts, yet require further evaluation using genomic data. Others (GMYC, 546 PTP) seemingly delimit populations, a conclusion in line with prior evaluations of simulated and empirical datasets. The latter was hypothesized to be impacted by inflated levels of inter-548 population divergence, a phenomenon which has been observed in other co-occuring desert fauna subjected to the same historical climatic fluctuation. These observations prompt us to 550 support the elevation of A. marmoratus, but additional groups withinthe A. tigris complex will requiring additional data to properly adjudicate.  Tables   786   Table 1   TAX=qualitative hypothesis derived from current taxonomy; CD5=CD-Hit with 5% divergence threshold for clustering; CD2= CD-834

Figures and
Hit with 2% divergence threshold for clustering; GMYC=Generalized mixed Yule coalescent; PTP=Poisson tree process; ROS=Clade delimitation based on statistical tests of monophyly. Node labels represent posterior probabilities from Bayesian inference (B), and 836 neighbor-joining bootstrap support (A; omitted for shallow nodes). Nodes passing the statistical test of reciprocal monophyly are indicated with a red dot. Clades in both trees are colored according to the ROS model. 838 Note that "tigris" refers to both A. t. tigris and A. t. punctilinealis, as these subspecies are not differentiated in the models. A. t. mundus was excluded due to the geographically restricted sampling presented herein. 858 860 Figure 7: Map depicting results of an 'Alleles in Space' analysis for members of the Aspidoscelis tigris complex. Relative genetic differentiation is interpolated across the landscape, with 862 individuals (as dots) colored by the CD5 model (see Fig. 1B). Heat map colors depict relative genetic distances per unit of geographic distance, with red corresponding to a greater genetic 864 distance, and blue being lesser genetic distance. 866