World Heritage lizard: population genetics and species status of the range-restricted Hamelin skink, Ctenotus zastictus

The Shark Bay World Heritage region in western Australia is home to a number of species of substantial conservation concern. Among these is a small scincid lizard, Ctenotus zastictus , which represents one of the most geographically-restricted vertebrates on the Australian mainland. The long-term persistence of Ctenotus zastictus is threatened due to the small size of its range, isolation from suitable habitat patches elsewhere, and potential impacts from climate change and mining. Accordingly, conservation efforts in Australia have targeted C. zastictus as the focus of protection. But this attention might be unwarranted – the species might not be evolutionarily unique. Previous genetic assessments have suggested limited differentiation between C. zastictus and its putative sister taxon, and the taxonomic status of C. zastictus has never been formally evaluated. Here, we use population genomic, phylogenetic, and ecoclimatic analyses to characterize the species status of C. zastictus in context of its closely-related congeners. In doing so, we explore the practical and conceptual challenges of revising species boundaries in threatened species, many of which are also rare and range-restricted. We demonstrate that C. zastictus is a coherent evolutionary unit that has been isolated from its putative sister species for at least two million years. Based on these results, we recommend that C. zastictus should retain its taxonomic status.


INTRODUCTION
Describing and revising species boundaries is not just a theoretical exercise (Avise 1989;May 1990;Hey et al. Devitt et al. 2019;Stanton et al. 2019). While some government and conservation agencies have moved towards more nuanced measures of biodiversity (e.g., evolutionary significant units; Moritz 1994;Pennock and Dimmick 1997), many still allocate resources according to assessments done using standardized taxonomies (Mace 2004). Thus, elevating the taxonomic status of an at-risk set of populations -to subspecies or species level -can immediately result in additional research attention, funding, and legal protections that may facilitate the long-term persistence of those populations (Morrison et al. 2009). Conversely, downgrading a population's taxonomic status can result in decreased resources and protection and hasten populations' decline and extirpation.
Thus, modifying taxonomy requires extra caution and care when the species involved are endangered or threatened (Hedin 2015). Some have argued that biologists should be conservative whenever changes to taxonomic status might put populations at increased risk of extinction (McCormack and Maley 2015). Others have countered that biologists should not retain taxonomically-dubious species (Zink et al. 2016). Doing so can directresources awayfrom other pressing priorities, including other species that might be under threat (Isaac et al. 2004), and can create apparent credibility gaps wherein science appears driven by underlying agendas rather than a consistent application of taxonomic standards (Zink et al. 2016). The middle ground might be to rigorously evaluate the case for and against changes to a species status. Unfortunately, there is no universally-accepted procedure for the robust delimitation of species (Hillis 2020;Burbrink and Ruane 2021;Hillis et al. 2021), in part because biologists do not agree on how to conceptualize species themselves (Stankowski and Ravinet 2021). Best practices include integrative approaches that quantify multiple axes of differentiation and the extent of reproductive barriers (Padial et al. 2010;Fujita et al. 2012;Pante et al. 2015b;Johnson et al. 2018), comparative approaches that place levels of differentiation in the focal taxa in context of patterns across close relatives (Tobias et al. 2010;Galtier 2019), and synthetic approaches that compare across multiple methods of species delimitation (Carstens et al. 2013). Ultimately, however, revised species boundaries -even when formalized into taxonomic revisions -remain a hypothesis subject to change with additional data (Hey et al. 2003;Isaac et al. 2004;Hedin 2015).
Here, we explore the challenges of species delimitation in threatened and endangered species by revisiting the species status of the vulnerable Hamelin skink, Ctenotus zastictus. Ctenotus, with more than 100 species, is perhaps the most diverse vertebrate genus in Australia (Cogger 2014). Many Ctenotus species are characterized by moderate to large geographic ranges (median size: ~140,000 km 2 , Prates et al. 2022) and are ecologically dominant in the communities in which they occur (Pianka 1969a(Pianka , 1969bRabosky et al. 2014a). However, C. zastictus is restricted to an extremely small region (~25 km 2 ; see Methods) found within the Shark Bay World Heritage Region in Western Australia (Fig. 1). This small range makes the species one of the world's most range-restricted non-insular lizards (Meiri et al. 2018). The species appears to be limited to only certain habitats; it is exclusively found on sandplains populated by spinifex grasses (Triodia sp.) and mallee Eucalyptus (Storr 1984). The taxon is absent from adjoining Acacia shrubland habitat (Cabrelli and Hughes 2015), and extensive surveys on adjacent sandplains have yet to document additional localities (McKenzie et June 2022al. 2000Aplin et al. 2008). Thus, the species range effectively consists of its type locality, although the species appears to be relatively abundant within that narrow region.
As might be expected given its small geographic range and specialized habitat, C. zastictus is at risk of extinction (Senior et al. 2021). The Australian Government currently lists C. zastictus as Vulnerable, and the International Union for the Conservation of Nature lists it as Near Threatened on the Red List of Threatened Species. This threat of extinction is exacerbated by climate change and potential mining impacts. Under current rates of climate change, C. zastictus is predicted to lose 99% of its habitat range by 2050 (Cabrelli and Hughes 2015). As the species range is only ~25 km², such a contraction would mean its range would completely disappear. In addition, C. zastictus partially occurs on Coburn station, which is being developed as a mining site. Because of its vulnerable status and these compounding effects, C. zastictus has been the center of conservation efforts. In 2015, the non-profit Bush Heritage Australia bought and restored a former pastoral station on which C. zastictus can be found, in part to conserve this species (Weber and Lewis 2016).
However, the taxonomic status of C. zastictus remains uncertain, and this focused attention might be unwarranted. The species was originally described based solely on morphological data; it purportedly has a unique pattern of spots relative to other species in the genus (Storr 1984). Although the vast majority of Ctenotus taxa were described based on body patterning (Cogger 2014), morphology is not always an accurate metric for distinctiveness in Ctenotus (Rabosky et al. 2014b. Across the genus, body markings are shared across species, and many species have polytypic patternings. Thus, genetic data have revealed that morphology and species identity are often mismatched, leading to both provisional and formal revisions of species boundaries across Ctenotus (Hutchinson et al. 2006;Kay and Keogh 2012;Rabosky et al. 2014bRabosky et al. , 2017Singhal et al. 2018b;Prates et al. 2022). Thus, the morphological distinctiveness of C. zastictus might not translate to Figure 1: Ctenotus zastictus shown with all known sampling sites from which vouchered specimens are known (Table S1). Photograph shows C. zastictus habitat (credit: DLR). C. zastictus is only known from very few unique localities. June 2022 evolutionary distinctiveness. Further, recent phylogenomic analyses have suggested that C. zastictus might be closely allied to-and potentially conspecific with-C. pallasotus. Ctenotus pallasotus was recently split from C. duricola ) and the species ranges across the Northwest Cape and western Pilbara regions of western Australia (Singhal et al. 2018b;Prates et al. 2022). Both C. pallasotus and C. zastictus share broadly similar patterns of lateral stripes and spots (Storr 1984;Rabosky et al. 2017). If C. pallasotus Doughty and Rabosky 2017 and C. zastictus Storr 1984 are conspecific, then C. pallasotus is a junior synonym of C. zasticus. As a result, the south Shark Bay populations that currently comprise C. zastictus would no longer have any special status at the species level, which might affect how its populations are treated in management and conservation plans.
In this study, we evaluate the evidence for the evolutionary distinctiveness of C. zastictus. We analyze genomic data for C. zastictus and its close relatives and use phylogenetic, population genetic, and ecoclimatic methods to determine its species status and evolutionary history.

Sample acquisition
We obtained full-precision georeferenced locality data for all vouchered C. zastictus specimens that could be located (n = 10) comprising seven from Western Australian Museum (WAM), one from the South Australian Museum (SAM), and two from the University of Michigan Museum of Zoology (UMMZ). Specimens held by the UMMZ were collected in 2013 under permit SF009161 from the Western Australia Department of Environment and Conservation, where lizards were captured using a combination of pitfall traps with drift fences.

Data collection and processing
We characterized patterns of genomic divergence across C. zastictus and its close relatives. Ctenotus is a species-rich clade of 107 recognized species (Uetz et al. 2021). The genus has traditionally been divided into a number of "species groups", although these have generally not been assessed from a phylogenetic perspective. Ctenotus zastictus is assigned to the atlas group, which spans 15 nominal species (Storr et al. 1999). We sampled all of these species: C. alacer, C. ariadne, C. atlas, C. decaneurus, C. duricola, C. dux, C. iapetus, C. impar, C. pallasotus, C. piankai, C. quattordecimlineatus, C. rhabdotus, C. xenopleura, C. yampiensis, and C. zastictus. To this, we added four species that our previous phylogenetic analyses suggested also fall within the atlas group: C. grandis, C. hanloni, C. maryani, and C. serventyi (Rabosky et al. 2014a;Singhal et al. 2017;Prates et al. 2022). For these species, we analyzed previously-published double-digest restriction-site associated DNA (ddRAD; n = 124; Table S2; Singhal et al. 2017;Prates et al. 2022) and mitochondrial DNA (mtDNA; n = 317; Table S2; Rabosky et al. 2014a). Additionally, we included ddRAD and mtDNA data from L. desertorum for use as outgroups.
To collect ddRAD data, we digested extracted DNA with the restriction enzymes PstI and MseI, barcoded fragments, and then selected for fragments between 150 -250 base pairs (Peterson et al. 2012). ddRAD libraries were pooled in sets of 100 and then sequenced using paired-end reads on either an Illumina HiSeq 2500 or 4000. We then generated two different assemblies and called genotypes using iPyrad v0.9.81 (Eaton and Overcast 2020). First, we were interested in evolutionary relationships among the atlas June 2022 group, and we thus assembled data across all sampled individuals (n = 124). Second, for the majority of our analyses, we focused on C. zastictus and the five species to which it is most closely-related: C. duricola, C. pallasotus, C. piankai, C. rhabdotus, and C. serventyi (n = 40; Fig. S1). Across both assemblies, we kept assembly and genotype calling parameters the same: clustering identity = 85%, minimum depth for base calling = 6, maximum percentage variant sites per locus = 20%.
To collect mtDNA data, we used Sanger sequencing to sequence the cytochrome B (cytb) gene, using primers from (Rabosky et al. 2009). Sequence chromatograms were assembled and then checked by eye using Geneious v2021 (Kearse et al. 2012). Sequences were aligned using mafft v7.471 (Katoh et al. 2009).

Data analysis
To determine the evolutionary distinctiveness and history of C. zastictus, we conducted a series of phylogenetic, population genetic, and ecoclimatic analyses.

Phylogenetic Inference
Because the atlas group was not defined based on a formal phylogenetic hypothesis, we first determined phylogenetic relationships among the species hypothesized to belong to the group. To do so, we generated a concatenated ddRAD alignment across all loci that were sampled for >60% individuals. We then used IQTree v2.1.3 (Minh et al. 2020) to infer an individual-level phylogeny and to calculate node support via approximate Shimodaira-Hasegawa likelihood ratio tests (SH-aLRT; Guindon et al. 2010). We also used IQTree to infer a mitochondrial gene tree, partitioning sites in cytb by coding position.
Based on these results (Fig. S1), we conducted a series of analyses focusing on C. zastictus and its five closest relatives: C. duricola, C. pallasotus, C. piankai, C. rhabdotus, and C. serventyi. We determined phylogenetic relationships among sampled individuals and taxa using three analyses. First, we generated a concatenated ddRAD alignment across all loci that were sampled for >50% individuals. Using this 9.9K locus and 1.9 Mb sites alignment, we inferred an individual-level phylogeny and 1000 ultrafast bootstraps with IQTree. In addition, we calculated nodal support as SH-aLRT values. Second, we used SNAPP v1.5.1 as implemented in BEAST 2 to infer a species tree (Bryant et al. 2012;Bouckaert et al. 2014). We randomly selected one single nucleotide polymorphism (SNP) per locus and only retained SNPs sampled in >80% of individuals. To reduce run times and improve convergence of SNAPP, we randomly subsampled three individuals per nominal species. We ran two parallel SNAPP analyses for 200K generations each (20% burn-in) and summarized results using TreeAnnotator. We then repeated this analysis across two additional subsampled datasets to ensure our results were robust to subsampling. Finally, we inferred a dated mtDNA phylogeny. We partitioned the cytb alignment by codon position and used ModelFinder as implemented in IQTree to determine the partitioning scheme with the highest likelihood and the best-fitting substitution model per partition (Lanfear et al. 2012). We used BEAST v2.6.4 for phylogenetic inference under a strict clock. We set a uniform prior on the substitution rate (0.7 x 10 -2 to 1.3 x 10 -2 substitutions per site per million years); this prior is informed by estimates of mitochondrial rate evolution across a diverse set of lizards and snakes from (Allio et al. 2017 (Skinner et al. 2011). We ran the MCMC analysis for 100e6 steps (20% burn-in), sampling every 5,000 steps, and we used Tracer to assess convergence and TreeAnnotator to summarize the run.

Population Clustering
We used two forms of population clustering analyses to determine how distinct individuals are in genotypic space. For all these analyses, we focused on C. zastictus and its five closest relatives (n = 37; three individuals were dropped due to missing data), and we used a dataset that randomly sampled one SNP from those loci that were >70% complete (n = 3822 SNPs). First, we ordinated all individuals using a principal component analysis (PCA; adegenet v2.1.4, Jombart 2008); this model-free approach is an effective way of seeing how distinct individuals are in genotypic space. Second, we conducted statistical population clustering using ADMIXTURE v1.3, which assigns individuals to genotypic clusters by minimizing levels of Hardy-Weinberg and linkage disequilibrium within a cluster (Alexander et al. 2009). For the ADMIXTURE analyses, we excluded C. serventyi and C. rhabdotus, because phylogenetic and PCA results suggested these two species were deeply divergent (Fig. 2,  3). We modeled the remaining data across 1 through 5 population clusters and measured model fit by the cross-validation error, which reflects how often individuals were misassigned to the wrong cluster. Further, because previous studies have shown that sampling imbalances across clusters can lead to biases in clustering (Puechmaille 2016), we randomly subsampled three individuals from each nominal species and then re-ran ADMIXTURE. We repeated this analysis an additional two times to ensure results were robust to subsampling.

Ecoclimatic Divergence
The analyses described above suggested that populations currently assigned to C. zastictus are most closely affiliated with those of C. pallasotus. We therefore conducted additional analyses to determine whether C. zastictus occurs in a distinct ecoclimatic space from its putative sister taxon. First, we compared where the two species occur in soil and climate space relative to each other and to their general geographic region. We randomly sampled 1000 points across a rectangle that bounds the geographic ranges of both species. For these random points, and the sampling sites for both focal species, we then extracted the local soil profile across 163 variables measuring properties such as silt, soil depth, soil pH, and water capacity (Grundy et al. 2015) and the local climate profile across 19 WorldClim variables at 5 arc-minute resolution measuring temperature and precipitation (Fick and Hijmans 2017). For the soil and climate datasets, we removed any variables that were missing across >25% of sampled sites and then ordinated sites using a principal component analysis (PCA; as implemented in the R package missMDA; Josse and Husson 2016). Second, we constructed an ecological niche model for C. pallosotus. Ctenotus zastictus has too few sampled points to allow accurate inference of a niche model. We used all known localities for C. pallosotus (n = 32; Fig. 2). We modeled the bioclimatic data only using MaxEnt v3.4 (Phillips and Dudík 2008), first fitting the model to all climatic variables and then re-fitting the model to only include the top variables (e.g., variables with permutation importance scores >5%). Our final model was built on six variables. We ran MaxEnt June 2022 with the default parameters, including a maximum of 500 iterations and 0.00001 convergence threshold. Adequacy of the model was evaluated both by review of the model's receiver operating characteristic (ROC) curve and predicted suitability under the model at known C. pallosotus localities. We extracted the niche suitability of C. zastictus localities under the C. pallosotus model as a measure of niche divergence.

Genetic assessment of species boundaries
We used two approaches to test whether patterns of genomic variation across the group are consistent with independently-evolving lineages. First, we explored a benchmark approach, in which the divergence of the putative taxa is compared to other recognized species in the clade (Tobias et al. 2010;Galtier 2019;Leaché et al. 2021). Using data from all atlas group species, we calculated nuclear d xy , nuclear Fst, and mitochondrial d xy across all pairwise comparisons of species. Per species and per metric, we then retained the lowest value as a species-specific estimate of divergence. Nuclear d xy is an absolute measure of genetic diversity that encompasses variation that is both segregating within populations and between populations. In contrast, Fst is a relative metric that is sensitive to the amount of variation within populations; reduced levels of within-population diversity can lead to increased Fst estimates (Cruickshank and Hahn 2014). Second, we used the statistical species delimitation approach implemented in *BFD (Leaché et al. 2014). Here, we tested two models: one in which C. zastictus was its own taxon and one in which it was lumped with C. pallasotus. Using the same SNP dataset used for SNAPP, we conducted a path-sampling analysis (100 steps of 200K length each) to calculate the marginal like-lihoods of these two models. We defined lambda broadly using a hyperprior (gamma distribution, (2, 200)) and defined the theta prior gamma distribution by setting alpha / beta equal to average sequence divergence among individuals within C. pallasotus. We then used Bayes Factors to compare likelihoods across models, after checking for convergence in Tracer. As we did with our SNAPP analysis, we repeated this analysis across two additional subsampled datasets to ensure our results were robust across subsamples. We recognize the assumptions and limitations of statistical species delimitation (Sukumaran and Knowles 2017;Leaché et al. 2019) and thus see this analysis as additional -but not conclusive -evidence regarding the species status of C. zastictus.

Demographic History
Our phylogenetic and population clustering results suggest that C. zastictus is evolving independently from its close relatives (see Results). Thus, we inferred the demographic history of these distinct populations and their close relatives. For all individuals sampled for C. zastictus and its sister species C. pallasotus, we randomly subsampled one SNP from each locus that was >70% complete (n = 3822 SNPs) and used this dataset to construct a two-dimensional folded joining site-frequency spectrum (2D-JSFS). We projected the 2D-JSFS down to 4 haplotypes for C. zastictus and 16 for C. pallasotus to account for missing data. We then used dadi v2.1.1 to fit a series of six divergence models to these two taxa (Fig. S2, Gutenkunst et al. 2009). These demographic models were all variants on a basic isolation-with-migration model. Previous work in Ctenotus has suggested that historical population expansion was likely prevalent across the genus (Singhal et al. 2017(Singhal et al. , 2021, which accords with expansion of arid June 2022 habitats in Australia during the late Miocene -mid-Pliocene (Pepper and Keogh 2021). Thus, we focused on demographic models that either modeled no population change or simultaneous population expansions only. Specifically, across our six models, diverging lineages experienced either no, symmetric, or asymmetric continuous migration, and in a subset of models, lineages experienced instantaneous population size change in the past.
To identify the best fitting model and its parameters, we followed the approach of Portik et al. (2017), running multiple rounds of optimizations multiple times. Because our SNP dataset consisted of unlinked SNPs, we calculated likelihoods using a standard approach and then determined the best-fitting model by calculating Akaike weights. We converted from coalescent units to real-time units using a published nuclear genome mutation rate for squamates (7.7 x 10 -10 ; Perry et al. 2018) and a generation time of 1 generation per year. Finally, we estimated parameter uncertainty by constructing 100 bootstrapped datasets and repeating the fitting approach for the best-fitting model only. Population-level sampling showed that C. zastictus forms a well-supported monophyletic group (Fig. 2, Fig. 3), with no evidence for nuclear-mitochondrial discordance (Fig. 3). Using a dated mtDNA phylogeny, we estimate that C. zastictus and its sister species C. pallasotus split ~3.8 million years ago (Fig. 2). A PCA of genotypic data reveals that C. zastictus formed a distinct cluster from C. pallasotus (Fig. 3B, S3). When these data were then used for statistical genetic clustering, the best model lumped C. zastictus and C. pallasotus into the same genetic cluster (Fig. 4, Table S3). However, sampling effort was uneven across the species analyzed (C. zastictus n = 3, all other species n = 8 -12) Given that sampling bias can affect genetic clustering (Puechmaille 2016), we randomly subsampled three individuals from each species three times. Across all subsampled datasets, the best-model identified each recognized species as its own genetic cluster (Fig. 4, Table S3).

Phylogenetic analyses of both the
We found that C. zastictus and its sister species C. pallasotus occupied adjacent   Table S3) and (B) subsampled dataset of three individuals per recognized species (best K = 4; June 2022 but unique regions in both soil and climatic space (Fig. 5A & 5B). Because we could not construct an ecological niche model for C. zastictus, we could only evaluate if C. zastictus localities are suitable for C. pallasotus. Our ecological niche model for C. pallasotus showed high suitability for known C. pallasotus localities, as expected (average suitability: 0.75; 95% range: 0.17 -0.99; Fig.  5C). However, C. pallasotus's niche model had extremely low suitability across C. zastictus's localities (average suitability = 0.009; range: 0.009 -0.009).
We considered two approaches to assessing species boundaries. First, we used a benchmark approach, where genetic divergence for C. zastictus was compared to those for recognized species across the atlas group. These results showed that C. zastictus had the lowest levels of absolute genetic divergence (as measured by d xy at the nuclear and mitochondrial genome) across all recognized species but its relative genetic divergence (as measured as Fst for the nuclear genome) was about average among recognized species (Fig. 6). Second, using a statistical approach to species delimitation, we found much greater support for a model in which C. zastictus and C. pallasotus were split as separate species rather than lumped together (Table S4). These results were consistent across three subsampled datasets.
Given evidence that C. zastictus is a distinct evolutionary unit, we inferred the demographic history between it and its sister species C. pallasotus. In the best-fitting demographic model, there were low levels of asymmetric gene flow after the two lineages split (Fig. S4, Fig. S5, Table S5). Divergence time between the two lineages was estimated at ~2.4 mya, and the C. pallasotus population was estimated to be ~12x larger than that of C. zastictus.

Species status of Ctenotus zastictus Storr 1984
We focused our phylogenetic, population genetic, and ecoclimatic analyses on determining whether Ctenotus zastictus shows evidence of being evolutionary cohesive and distinct from its close relatives in the scincid genus Ctenotus. Our results indicate that C. zastictus comprises a coherent evolutionary species, with a substantial history of isolation from its presumed sister taxon. Accordingly, we argue that the current taxonomic status of C. zastictus should be retained. Our phylogenetic and population genetic analyses confirm that C. zastictus is a distinct and independent lineage. Ctenotus zastictus forms a highly-supported monophyletic clade in phylogenies constructed with both nuclear and mitochondrial data (Fig. 2, Fig. 3). In the species description for C. zastictus, Storr hypothesized that the species was sister to C. iapetus based on shared morphologies (Storr 1984). Instead, both mitochondrial gene trees and phylogenomic phylogenies reveal that C. zastictus is sister to C. pallasotus, a species that was recently split from C. duricola based on patterns of genomic and morphological distinctiveness ). The population genetic data mirror the phylogenetic data. In genotypic space, C. zastictus forms a cohesive group that is distinct from other closely-related species, including its sister species C. pallasotus (Fig. 3B, S3). When the same genotypic data are analyzed using genetic clustering approaches, C. zastictus is inferred to be its own unique genetic population, provided that sampling is even across lineages (Fig. 4, Table S3). Further, C. zastictus appears to be evolving independently. Across both phylogenetic and population genetic data, we found no evidence for admixture or genealogical discordance, suggesting either no or limited introgression between C. zastictus and its closely-related species.
There is also evidence that Ctenotus zastictus is ecologically distinct. It occurs in unique ecoclimatic space from its sister species C. pallasotus ( Fig. 5A-B), and niche modeling suggests the two species likely also have divergent climatic niches (Fig.  5C). Further, the two species are found in two different bioregions within Australia. Ctenotus zastictus is found in the Wooramel subregion of the Carnavon Basin, and C. pallasotus occurs in the Pilbara region (Thackway and Cresswell 1995). These regions circumscribe habitats within Australia and accurately predict occurrence and range limits of species, reinforcing that these two species live in distinct biotic gence seen for a given species. C. zastictus is shown in red. C. zastictus shows the lowest values of absolute sequence divergence (d xy ) among the recognized species in this group, though its relative sequence divergence (Fst) is on par with other recognized species. June 2022 and physiographic regions. However, the small range size and limited sampling for C. zastictus prevents us from quantifying exactly how distinct its niche is from its sister species. Both informal and formal approaches to species delimitation agree with these more descriptive analyses of C. zastictus cohesiveness and distinctiveness. If we use a benchmark approach, C. zastictus is similarly distinct as other recognized species in the atlas clade across relative measures of genomic divergence (Fst; Fig. 6). Across absolute measures (d xy at the nuclear and mitochondrial genome), C. zastictus is far less divergent than other species. The small population size of C. zastictus is the likely source of the discrepancy between patterns of relative and absolute sequence divergence. Ctenotus zastictus has low levels of genetic diversity segregating within its populations (Fig. S6), which then leads to higher levels of relative genetic divergence, whether measured by Fst or d a (Cruickshank and Hahn 2014). Given our discrepant results across relative versus absolute measures of sequence divergence, our benchmark approach is somewhat equivocal about C. zastictus's species status. In contrast, our statistical approach to species delimitation was decisive. Using the species delimitation approach *BFD, modeling C. zastictus as a distinct species was much more supported than a model in which it was lumped with C. pallasotus (Table S4). However, statistical species delimitation approaches can delimit populations as species (Sukumaran and Knowles 2017;Singhal et al. 2018a;Leaché et al. 2019) and thus these results should be interpreted cautiously.
Ctenotus zastictus was delimited based solely on morphological grounds, which have been shown to be unreliable in the genus Ctenotus (Rabosky et al. 2014bSinghal et al. 2017;Prates et al. 2022). By using genetic and ecoclimatic analyses and placing the evolution of C. zastictus in context of the broader clade, we confirm this initial morphological diagnosis and find that the species is evolutionary distinct and cohesive. Our results suggest that populations assigned to C. zastictus (south Shark Bay) comprise a distinct gene pool with a considerable (2+ million year) history of isolation from other Ctenotus populations, most notably the northern populations currently assigned to C. pallasotus. Importantly, we did not measure additional aspects of species ecology, morphology, or the extent of reproductive barriers, any of which could change our understanding of how distinct these lineages are. Accordingly, this conclusion -like all conclusions about species status -remains provisional (Fujita et al. 2012). That said, based on existing data, we recover no evidence that C. zastictus is taxonomically dubious. While it is less divergent than many other species in the clade (Fig. 6), it still shows strong evidence of being on an independent evolutionary trajectory. Thus, C. zastictus satisfies the criteria for recognition under an evolutionary species concept, and it should be continued to be recognized as taxonomically distinct from C. pallasotus and other Ctenotus species.

The origins of Ctenotus zastictus
Based on our demographic modeling, C. zastictus and its sister species C. pallasotus likely split nearly 2.4 million years ago, while estimates based on our dated mitochondrial tree suggest that the divergence time for the mtDNA lineages of these two species was 3.8 myr. Across both dating exercises, we neither have fossil data nor do we have an estimate of the mutation rate from within the genus. Instead, we rely on secondary calibrations (Skinner et al. 2011;Rabosky et June 2022al. 2014a) and estimates of mutation rate from other squamates (Allio et al. 2017;Perry et al. 2018) in order to convert our coalescent estimates of divergence time to actual time units. Thus, these dates should be treated as provisional, although they strongly suggest a pre-Pleistocene or early Pleistocene split between C. zastictus and C. pallasotus. The best demographic model further suggests that the two lineages exchanged migrants at low levels through divergence (0.02 -0.13 Migrants / generation; Table S5, Fig. S5). Historical introgression is common through species divergences in vertebrates (Pinho and Hey 2010;Mallet et al. 2016), and even low levels of introgression can affect divergence if migration is prolonged, populations exchanging migrants are small, or adaptive variants are shared between diverging populations. However, we found no evidence that this limited introgression is eroding divergence between C. zastictus and C. pallasotus.
The Shark Bay region and adjacent sandplains are home to a number of endemic taxa and/or taxa with disjunct or highly-differentiated populations (Storr and Harold 1978;Hopper and Gioia 2004;Edwards et al. 2007;He et al. 2013). Some species have distinct ranges across the western Australian coast between the Northwest Cape and Geraldton, whereas many other species have highly-restricted, predominantly coastal distributions in this region (Rabosky et al. 2004;Maryan et al. 2015). Cryptic population differentiation has been observed over short geographic differences along Shark Bay even for relatively high-vagility avian taxa (Walsh et al. 2021). These and other biogeographic patterns might be due to Pleistocene sea level changes: lower sea levels during glacial maxima would have exposed vast areas of coastal habitat -much of it sandplainand potentially re-connecting populations that are currently isolated by surficial geomorphology (Kendrick et al. 1991). In particular, much or all of today's shallow and saline Shark Bay marine environment would have been dry land during periods of reduced sea level, providing a possible dispersal corridor between present-day distributions of C. zastictus and C. pallasotus. However, our demographic modeling implies that expanded coastal habitat during the Pleistocene resulted in little gene flow between these species. Our results are consistent with increasing evidence that pre-Pleistocene physiographic processes have played a more important role than Pleistocene processes in mediating population divergence and speciation along the western coastal margin of Australia (Edwards 2007).

Species boundaries in rare and restricted taxa
Determining species status is notoriously difficult (Carstens et al. 2013;Pante et al. 2015a;Cadena and Zapata 2021), and it is even more difficult when the populations under consideration are rare or range-restricted. Further, many rare or range-restricted species-like C. zastictus-are vulnerable and endangered, which makes species delimitation a "high-stakes" affair (Hedin 2015;Devitt et al. 2019). Most species delimitation methods rely on distinguishing between infraspecific and interspecific patterns of variation (Lim et al. 2012). However, characterizing within-species patterns of genetic and phenotypic variation requires broad sampling, which is often challenging in rare or range-restricted species. Practically, this means that researchers working on rare or range-restricted species cannot easily use many of the tools commonly employed in species delimitation studies. We faced these challenges with C. zastictus. For example, we could not construct an ecological niche June 2022 model for C. zastictus and formally compare divergence across niche models between C. zastictus and C. pallasotus (c.f., Warren et al. 2008). As another example: one way to infer the extent of reproductive barriers is by looking for evidence of introgression between lineages. We both looked for mito-nuclear discordance across phylogenies and admixture in population clustering results as evidence of introgression. However, these results are inherently limited because we only sampled three individuals from C. zastictus. Finally, when rare species are delimited relative to more widespread species, sampling across the species can be imbalanced, leading to additional issues. For example, our sampling imbalances affected the number and identity of the clusters recovered using genetic clustering algorithms (Fig. 4). However, our work also showcases the potential to delimit rare species in the genomic age. Because collecting thousands of loci from a single individual is now affordable (Peterson et al. 2012), using only a few individuals per species, we can robustly measure basic metrics of genetic diversity and divergence (Fig. 6, Fig.  S6) and apply coalescent approaches to infer evolutionary relationships among species, to identify species boundaries, and to estimate demographic history (Fig. 3C, Table S4, Fig.  S4).