Population genomics of flat‐tailed horned lizards (Phrynosoma mcallii) informs conservation and management across a fragmented Colorado Desert landscape

Phrynosoma mcallii (flat‐tailed horned lizards) is a species of conservation concern in the Colorado Desert of the United States and Mexico. We analysed ddRADseq data from 45 lizards to estimate population structure, infer phylogeny, identify migration barriers, map genetic diversity hotspots, and model demography. We identified the Colorado River as the main geographic feature contributing to population structure, with the populations west of this barrier further subdivided by the Salton Sea. Phylogenetic analysis confirms that northwestern populations are nested within southeastern populations. The best‐fit demographic model indicates Pleistocene divergence across the Colorado River, with significant bidirectional gene flow, and a severe Holocene population bottleneck. These patterns suggest that management strategies should focus on maintaining genetic diversity on both sides of the Colorado River and the Salton Sea. We recommend additional lands in the United States and Mexico that should be considered for similar conservation goals as those in the Rangewide Management Strategy. We also recommend periodic rangewide genomic sampling to monitor ongoing attrition of diversity, hybridization, and changing structure due to habitat fragmentation, climate change, and other long‐term impacts.

phylogenetic effects of range size and abundance (Lawton, 1993), the mechanisms underlying these patterns remain poorly understood.For the agencies tasked with proactively conserving habitats and biodiversity, molecular data can inform realistic and practical management strategies to prevent rare species from reaching the point of collapse.Here, we apply genomic data to the management of a rare species inhabiting the Colorado Desert of southwestern North America, a region containing extensive natural areas but also growing cities, agriculture, renewable energy developments, and military bases.

| Study system
Flat-tailed horned lizards (FTHLs), Phrynosoma mcallii (Hallowell, 1852), are endemic to the Colorado Desert, specifically the Salton Trough (Figure 1).This small range-the smallest of 10 Phrynosoma species in the United States-spans the states of California and Arizona in the United States, and Sonora and Baja California in Mexico (Sherbrooke, 2003).Phrynosoma mcallii is nested within the Anota clade of Phrynosoma, which includes P. solare of the Sonoran Desert uplands, and the P. coronatum complex of coastal California and Baja California (Leaché & Linkem, 2015;Leaché & McGuire, 2006).The Colorado Desert is ecologically known as the Lower Colorado River Valley subdivision of the Sonoran Desert (Dimmitt, 2015) but is often treated as a distinct desert, particularly in California, because it lacks the regular summer precipitation and columnar cacti of the relatively lush uplands of Sonora, Baja California, and Arizona (Pavlik, 2008).Geologically, the Salton Trough, more than 5400 km 2 of which is below sea level, is a tectonically active fault-bounded structure, part of the depression flooded by the Gulf of California (Stokes et al., 1997).The trough is traversed by both the San Andreas fault system and the lower reaches of the Colorado River, which together have added a unique dimension to the region's geology and hydrology (Blount & Lancaster, 1990;Stokes et al., 1997).Periodic highstands caused by westward diversion of the Colorado River into the Salton Trough, triggered by either large floods or tectonic activity, resulted in Pleistocene-Holocene Lake Cahuilla, with an estimated highstand surface area of 5700 km 2 (Waters, 1983).In 1905In -1906, Lake Cahuilla was reincarnated as the Salton Sea when floodwaters caused irrigation canals to overflow, filling the lakebed once again (Gobalet & Wake, 2000).
Being endemic to the Colorado Desert has a cascade of effects on the life history of P. mcallii that may predispose it to rarity.
Given its latitude, elevation near or below sea level, and the rainshadow effect, the Colorado Desert's climate is extremely hot and arid; air temperatures may exceed 49°C, with surface temperatures approaching 71°C, and average annual rainfall in the driest sites is less than 76 mm (Dimmitt, 2015).Increasing global temperatures may pose a direct threat to this species, as daytime activity periods are decreasing (Sinervo et al., 2010).Ants, especially Messor and Pogonomyrmex, make up >97% of the diet, which is more than in other Phrynosoma species (Pianka & Parker, 1975;Turner & Medica, 1982).Prey availability is correlated with seed production by annual plants in spring, which is impacted by the amount and timing of winter precipitation (Beatley, 1967;Brown et al., 1979).
Some P. mcallii populations may fluctuate with drought cycles (FTHL ICC, 2003), but studies at the northwestern range limit indicate that populations negatively correlate with the density of an invasive plant (Brassica tournefortii), which occurs in higher abundance during wet years (Barrows et al., 2021;Barrows & Allen, 2009).To find adequate prey in the open desert, lizards need relatively large individual home ranges, up to 35,000 m 2 (FTHL ICC, 2003).More movement is needed for foraging and mating in larger home ranges, exposing lizards to mortality due to predation, the intensity of which also varies annually (Young & Young, 2000).FTHLs also produce relatively small clutches, compared to other Phrynosoma species, of three to seven eggs (Howard, 1974;Pianka & Parker, 1975).The extreme temperatures and aridity in the Colorado Desert are hypothesized to contribute to their low fecundity compared to other Phrynosoma species (Howard, 1974).
Several other aspects of the ecology and behaviour of FTHLs increase their sensitivity to human disturbance.Like other horned lizards, they rely on crypsis more than speed to evade predators and often freeze in response to danger, making them vulnerable to vehicles (FTHL ICC, 2003).Predation has been shown to be a top driver of population numbers (Goode & Parker, 2015;Muth & Fisher, 1992;Young & Young, 2000).Power lines and palm savannah provide artificial perches for avian predators (FTHL ICC, 2003;Leavitt et al., 2015); the link between development and increased predation is also well established with common ravens (Corvus corax) and desert tortoises (Gopherus agassizii) in the nearby Mojave Desert (Kristan & Boarman, 2003).As of 2003, the species has lost approximately 43% of its original habitat in the United States (FTHL ICC, 2003).Other potential threats include widespread habitat loss and fragmentation due to agriculture, urban expansion, off-road vehicle use, renewable energy development, and military activities (Barrows & Allen, 2009;USFWS, 1993).
Beginning in the early 1990s, a diverse group of state and federal stakeholders began to form conservation strategies that would maintain current population levels and avoid further regulatory restrictions on their land usage.This effort resulted in the 1997 Conservation Agreement signed by these stakeholders to '… conserve [P.mcallii] by reducing threats to the species, stabilizing the species' populations, and maintaining its ecosystem'.Resulting from this significant voluntary effort, the signatory agencies and offices then created the flat-tailed horned lizard rangewide management strategy (FTHL RMS;FTHL ICC, 1997, 2003), intended to maintain extant populations of the species in perpetuity.The FTHL RMS established five Management Areas (MAs) on both sides of the Colorado River at the core of the range, and one Research Area (RA) in California, that together cover 2274 km 2 of habitat managed by state, federal, military, and private stakeholders (Figure 1).(Mulcahy et al., 2006).However, mtDNA is essentially a single locus and cannot account for gene flow caused by male dispersal; a genome-wide multi-locus approach is preferable to robustly test demographic hypotheses in a coalescent framework (Gottscho et al., 2017).Here, we use a population genomics approach to further inform management strategy on state and federal stakeholder lands by addressing the following questions: 1. What is the population structure of P. mcallii and how are populations related phylogenetically?
2. Are there any disparities in genetic diversity across the landscape?
3. Are there significant barriers to gene flow among populations?
4. What are effective population sizes, and is there evidence of population size changes?
From a strategic perspective, determining how closely the current management units of the FTHL RMS correspond to biological populations has direct implications for allocating resources and prioritizing conservation and management actions for respective populations and habitats.

| Genomic DNA samples
We extracted genomic DNA from tissue samples of P. mcallii collected throughout its known range using the PureLink Genomic DNA kit (ThermoFisher Scientific, Waltham, MA).Most of these tissue samples were in the form of tail clips, toe clips, or blood samples that were collected by various individuals (see Acknowledgements) and were preserved in 96%-100% ethanol and stored at −20 or −80°C.We then screened samples for quality (high nucleotide concentration and molecular weight).In our final dataset, we included 45 high-quality samples (Table 1), 42 representing P. mcallii spanning the extent of its known range, one representing the outgroup species P. platyrhinos (the desert horned lizard), and two putative P. mcallii × P. goodei hybrids (see Mulcahy et al., 2006).

| Library preparation and sequencing
We collected DNA sequence data using the double-digest restriction-associated-DNA sequencing (ddRADseq) method (Peterson et al., 2012)

| Bioinformatics
The Python pipeline ipyrad v0.9.84 (Eaton & Overcast, 2020) was used for quality control, reference-based assembly of RAD loci, genotype calling, and generating alignments of loci and SNPs.The P. platyrhinos genome (1.9 Gb long, containing 6 macrochromosomes and 11 microchromosomes) was used as the reference (Koochekian et al., 2022).Samples were demultiplexed by their unique barcode and adapter sequences, implementing a strict filter for adapter and primer sequences.Sites with Phred quality scores under 33 (99.95%) were converted to 'N' characters, and reads with >5 N's were filtered out.
The filtered reads were clustered using a threshold of 90%.Additional filters applied to consensus sequences include a maximum of 50% heterozygous sites per locus, a maximum of 5% uncalled bases in consensus sequences, and a maximum of 5% heterozygotes in consensus sequences.Separate sequence and SNP matrices were assembled with (n = 45) and without (n = 42) the outgroup sample of P. platyrhinos and the two putative P. mcallii × goodei hybrids.To control the level of missing data, we set the minimum number of samples per locus for output to 21 for the ingroup-only and 24 for all samples.The ingroup matrices (n = 42) were subjected to additional filtering steps using VCFtools v0.1.13(Danecek et al., 2011) to better meet the assumptions of downstream population genetic and demographic modelling analyses.The minor allele frequency threshold was set to 0.05, and a maximum of two alleles were allowed per locus.A sliding window was used to thin variants such that no more than one variant was allowed per 10,000 bp window.The maximum missing data threshold was set to 50%.

| Phylogenetic network and tree
We constructed a genetic network with SplitsTree v4.16 (Hudson & Bryant, 2006) from the concatenated RAD loci for the ingroup samples (n = 42; 942,435 bp).Network methods are appropriate for population-level studies where non-bifurcating relationships and admixed samples are expected (Posada & Crandall, 2001).The network was constructed using uncorrected 'p' distances and the NeighborNet algorithm (Bryant & Moulton, 2004).We also constructed a maximum likelihood phylogeny with IQ-TREE v2.1.3(Nguyen et al., 2015) of the concatenated data that included the outgroup and hybrids (n = 45; 875,955 bp).We used ModelFinder (Kalyaanamoorthy et al., 2017) to test 286 DNA substitution models in IQ-TREE.The best-fit model (K3Pu + F + R2) was selected following the Bayesian information criterion (BIC), and 1000 bootstrap replicates were performed.

| Population structure
We estimated population structure using multiple approaches.In all cases, the ingroup dataset (n = 42) contained only filtered SNPs as described above.First, we used discriminant analysis of principal components (DAPC) implemented in the R package adegenet v2.1-4 to cluster individuals without any prior grouping (Jombart et al., 2010).We determined the optimal number of clusters (K) by running a principal components analysis (PCA), retaining all of the TA B L E 1 Lizards sampled in this study and their approximate sampling locations.PCs ( 41), and using the BIC to select the best fitting model (testing K values 2-10).We retained 18 PCs for the subsequent DAPC analysis based on the results of a-score optimization and spline interpolation (Jombart et al., 2010).
We also used the maximum likelihood method Admixture v1.3.0 (Alexander et al., 2009) to determine the optimal number of populations (K), population structure, and admixture proportions.The cross-validation error (CVE) was used to select the optimal number of populations (K), testing values ranging from 2 to 10.The analyses were repeated 10 times to measure uncertainty in CVE.We also used Admixture to calculate the fixation index (F ST ) between populations.Admixture bar plots were generated with Structurly 0.1.0(Criscuolo & Angelini, 2020).

| Estimated effective migration
We used the programme EEMS (estimated effective migration surfaces) to model the relationship between genetics and geography (Petkova et al., 2016).This approach uses a genetic distance matrix and a geographic distance matrix to ascertain deviations from a null isolation-by-distance model.The method estimates withinpopulation effective diversity (q) and among-population effective migration (m) parameters.Deviations from the isolation-by-distance model are visualized in the form of a spatial heat map displaying barriers and corridors for migration and local diversity hotspots.We calculated an Euclidean distance matrix based on the same ingroup SNP matrix mentioned above.We overlaid a grid of 300 demes in the shape of the approximate range of this species.We ran the programme for 20 million steps, discarding the first 2 million as burn-in, sampling every 25,000 steps.

| Demographic modelling
We tested demographic models for the divergence of P. mcallii across the Colorado River by simulating the two-dimensional (2D) joint site frequency spectrum (JSFS) of genetic variation between populations using the software package moments (Jouganous et al., 2017).
We optimized 16 2D demographic models (Figure S1) that parameterize isolation, migration, and effective population size (Leaché et al., 2019;Portik et al., 2017).We focus on demographic models with parameters that pertain to the allopatric and parapatric modes of divergence worth exploring in a terrestrial lizard distributed across a river barrier.The simplest model includes divergence without migration, while the most complex models include multiple epochs, asymmetrical migration, secondary contact, and instantaneous population size change.We compared models by dividing the samples into two populations located on opposite sides of the Colorado River, which correspond to northwestern and southeastern populations.
We projected down to smaller sample sizes to maximize the number of segregating sites used in analyses (northwest = southeast = 20) using the programme easySFS (Gutenkunst et al., 2009).Four rounds of model optimization were performed using the log_fmin optimizer (Jouganous et al., 2017).Each round of model optimization used 100 replicate searches, 50 iterations per step, and a multinomial approach to estimate the log-likelihood of each model.The parameters from the best-scoring replicate were used as starting values for the next round of optimization.After the final round of optimization, the replicate with the highest likelihood for each model was used to calculate AIC scores, ΔAIC scores, and Akaike weights (Burnham & Anderson, 2002).Replicate analyses were conducted to ensure that the optimization routine was stable.Demographic model fitting was performed using the moments pipeline (Leaché et al., 2019).
To convert the unscaled population parameters to demographic terms, we assumed a generalized lizard mutation rate (μ) of 5.17e −9 substitutions/site/generation (Bergeron et al., 2023) and a generation time of 1 year (Howard, 1974).For the time parameter T, we used the equation we used N ref = θ/4 μL, where μ is the mutation rate, L is the number of loci multiplied by their length, and θ was derived from the moments output.To produce 95% credible intervals for the estimates, we bootstrapped the top-ranked model (symmetric migration with population size change 'sym_mig_size').We generated 50 new bootstrapped replicate datasets, each containing the same number of SNPs and samples, using sampling with replacement of the unlinked SNPs.Each of the 50 bootstrapped datasets was then run in moments using the same procedure and settings as for the model selection procedure described above.We used R to post-process the bootstrap results to calculate summary statistics for each of the model parameters (mean and 95% credible interval).We also made density plots to compare the bootstrapped to the optimized values.

| Phylogenetic network and tree
The phylogenetic network of concatenated RAD loci for the ingroup

| Population structure
The BIC and CVE scores obtained from DAPC and Admixture indicate that a two-population model best fits the data (Figure S2).The K = 3 and K = 4 models, which are only marginally less supported than K = 2, provide further geographic structure (see below for detailed descriptions of the geographic structure under these two models).with mixed ancestry, the results from DAPC and Admixture were congruent in delimiting population structure.F ST values between populations calculated by Admixture were low-under the twopopulation model, F ST was 0.17, while under the three-population model, F ST was 0.18 between the eastern and southwestern populations, 0.37 between the eastern and Coachella populations, and 0.29 between the southwestern and Coachella populations.These F ST values indicate that the majority of genetic variation is shared between these populations.

| Estimated effective migration
The results of the EEMS analysis show that the lower Colorado River and the Salton Sea represent the major geographic barriers to migration among demes (Figure 5).

| Demographic modelling
Ranked models are shown in Figure S1 and Table S2.

| DISCUSS ION
In this study, we present analyses of novel genome-wide data of a rare desert lizard species with a restricted range, with the overarching goal of informing the rangewide management strategy (RMS).
In addition to corroborating published studies documenting divergence across the lower Colorado River in this species and others, we find strong evidence of a recent, severe bottleneck.Interpreting and applying our results requires consideration of other sources of evidence, including geology and hydrology, as the Colorado Desert is a region of environmental extremes, shifting from desert to freshwater environments during the Quaternary, and undergoing anthropogenic transformation in the past century.

| Phylogenetic and clustering analyses confirm a Colorado River break and northwestern range expansion
The first objective of this study was to determine the rangewide population structure and phylogeny of P. mcallii as a starting point to guide downstream analyses and management.The unrooted SplitsTree network (Figure 2) and rooted concatenated phylogeny

| Mapping disparities in genetic diversity across a fragmented Colorado Desert landscape
To model the relationship between genomics and geography, we used the EEMS method (Petkova et al., 2016), which is designed for a study with geo-referenced individuals and with population structure that is broadly consistent with isolation by distance.
EEMS produces maps of the effective diversity rate (q), the expected genetic dissimilarity of two individuals sampled from a given location.The most prominent hotspots of diversity were found to be west of the Salton Sea, south of Yuma, and at El Pinacate.Areas with lower diversity include the Coachella Valley and southernmost Colorado River delta in Sonora (Figure 5).
These results are consistent with mtDNA data, which also show that P. mcallii populations west of the Colorado River have less genetic differentiation compared to those southeast of the river (Mulcahy et al., 2006).The low diversity in the Coachella Valley could be a result of natural range expansion in prehistoric times or modern habitat fragmentation from rapid development of the valley over the last century.Again, our results are consistent with published mtDNA data showing significant population structure, limited gene flow, and isolation associated with the Colorado River (Mulcahy et al., 2006).

| Modelling migration rates, effective population sizes, and population size change
After testing 16 divergence models of varying complexity, with and without migration parameters, the best-fit model included bidirectional migration across the Colorado River, albeit at very low levels (Table 2).We present the first genome-derived estimates of effective population size, with obvious relevance for this species of special concern (Table 2).Strikingly, we detected severe population declines (>99% in both populations) that occurred very recently.The best estimate for the initial divergence of the two populations across the Colorado River is late Pleistocene (95% CI 13-345 kya), while the age of the bottleneck is entirely within the Holocene (95% CI 176-4990 years ago).
What could explain the timing of the initial divergence and bottlenecks in both populations?If we assume the inferred dates are correct-and they depend on the assumed mutation rate and have additional sources of uncertainty not captured by bootstrapping, including model selection and parameterization-then one hypothesis is that the divergence or bottleneck occurred due to the filling of Lake Cahuilla (a lacustrine interval).In this region of environmental extremes, Lake Cahuilla formed periodically in response to the western diversion of the Colorado River in the Salton Trough, triggered by large floods or tectonic events (Stokes et al., 1997;Waters, 1983).This is a unique mode of origin that distinguishes this system from other lakes in the region, for example, in the Mojave Desert, that instead responded to fluctuating precipitation and temperature during Pleistocene glacial cycles.In the last 2000 years, there were four highstands of Lake Cahuilla, the most recent ending in AD 1430 (Waters, 1983).The youngest palaeo-shoreline is still visible at 12 m elevation, implying an estimated surface area of 5700 km 2 and a maximum depth of 95 m.Older shorelines are documented well into the Pleistocene, although their shorelines are more fragmented and deformed by tectonics (Waters, 1983).By the time Europeans arrived, this instance of the lake had vanished, implying disappearance in about 60 years.However, there are two difficulties in attributing a direct role to any specific lacustrine interval to either the initial divergence or the bottleneck.First, the 95% CI for both events spans multiple lacustrine intervals.Second, the bottleneck event was equally severe on both sides of the Colorado River.Another hypothesis is that the bottleneck can be attributed to the dramatic anthropogenic habitat transformation that occurred in the past century.The 95% low estimate for the bottleneck is 176 years ago, barely earlier than large-scale habitat modification, so strictly interpreted, we would reject this hypothesis.However, it is likely that there are additional sources of uncertainty surrounding our divergence date estimates, so we believe that we cannot reject the hypothesis that the severe bottleneck resulted from the Anthropocene impacts described in the Introduction.

| Comparisons with co-distributed species
Comparative phylogeography can reveal underlying processes that influence shared distributional patterns.In the Colorado Desert, the most obvious physical barrier is the Colorado River, which is associated with stronger genetic differentiation in reptiles, arachnids, and mammals relative to flying insects and birds (reviewed by Dolby et al., 2019).Our results of low but statistically significant gene flow (1.9-30 migrants per generation) add to the evidence that Colorado River functions as a leaky filter barrier to migration for many nonvolant species-there is enough of a barrier to promote divergence, but post-divergence gene flow remains common (Dolby et al., 2019).A major challenge in directly comparing results from published studies is confounding factors from phylogeny, geographic sampling, choice of molecular markers, and model selection.For example, in comparisons across the Colorado River cited above, original studies used varying approaches to distinguish between the effects of gene flow and retained ancestral polymorphism, or to determine the directionality of gene flow (Dolby et al., 2019).For that reason, it is especially fruitful to compare the patterns found in this study to those of fringe-toed lizards (Uma), for which a highly comparable ddRADseq dataset is also available, including similar geographic sampling and near-identical library preparation, restriction, and size selection steps (Gottscho et al., 2017).Both P. mcallii and Uma are in the clade Phrynosomatidae (Leaché et al., 2015;Reeder & Wiens, 1996), and the two occupy similar habitats, with mostly overlapping ranges entirely within the Colorado Desert / Salton Trough (Gottscho et al., 2017).Both have a northwest-from-southeast nested phylogenetic pattern, divergence dates across the Colorado River in the Pliocene-Pleistocene, and detectable gene flow across the river.
Both have higher genetic diversity south and east of the Colorado River and lower diversity in the Coachella Valley.Given their relatedness and similar life histories, these similarities are not altogether surprising, but nevertheless collectively represent a critical validation of our results.
There are also key differences, however.Unlike P. mcallii, the Uma notata complex consists of at least four species across the  et al., 2020;Gottscho et al., 2017;Trépanier & Murphy, 2001).The distributions of these species closely correspond to the ranges of P. mcallii populations under K = 4 (compare Figure 3 in this study to figure 1 in Gottscho et al., 2017).This raises the question as to whether any identified P. mcallii populations have diverged enough to be considered (sub)species.With the caveat that this study was never intended to be a taxonomic revision, based on the current evidence, we think that they do not.As in this study, Gottscho et al. (2017) used BIC and CVE scores to determine optimal K values for downstream species delimitation.Excluding the species occupying the Mojave Desert (U. scoparia) and Mohawk Dunes (U. thurmanae), where P. mcallii is absent, 3-4 Uma clusters from Coachella to El Pinacate scored best (Gottscho et al., 2017, Table 2), whereas in this study, K = 2 was best supported and only marginally better than K = 1 (Figure S2).Uma inornata, endemic to the Coachella Valley, is also morphologically diagnosable and has been recognized as a distinct (sub)species well before molecular data (Heifetz, 1941).Our interpretation is that Uma shows more differentiation, structure, and ultimately speciation than P. mcallii across the same landscape because fringe-toed lizards are restricted to aeolian sand, for which they possess specialized adaptations, and which is patchily distributed.By contrast, P. mcallii occupies both sandy and intervening hardpan substrates (Rorabaugh & Young, 2009), increasing overall suitable habitat and gene flow between metapopulations relative to Uma-at least before modern habitat fragmentation.and State Ecological Reserve (Barrows et al., 2021;Barrows & Allen, 2009).By contrast, Uma inornata, a federally threatened endemic species, still occurs at multiple sites across the valley (Barrows & Heacox, 2021).Vandergast et al. (2016)  with one species or the other (Mulcahy et al., 2006).Additionally, analyses of nuclear DNA sequence data showed these same putative hybrids to be polymorphic at several sites in two genes (Leaché & McGuire, 2006).Furthermore, phylogenetic analyses of horned lizards based on mtDNA place P. platyrhinos (and P. goodei) sister to P. mcallii, whereas nuclear data place P. platyrhinos-P.goodei with P. modestum (Doliosaurus) and P. mcallii is placed with P. solare and the P. coronatum complex (Anota) (Leaché & McGuire, 2006).Here, genome-wide molecular data support the hybridization hypothesis between P. mcallii and P. goodei as the two hybrid individuals included were phylogenetically intermediate between these two species (Figure 3).

| Recommendations and next steps
The ultimate goal of this study was to use genomic data to inform the FTHL RMS, and in turn to sustain populations of P. mcallii in perpetuity.Here are four action items based on our findings that RMS stakeholders could consider for ongoing conservation efforts.Mulcahy et al., 2006), where the substrate transitions from pebble to sand and the two species narrowly overlap (Young & Young, 2000); additional hybrids with P. platyrhinos were reported near Ocotillo, CA (Stebbins, 2003).Natural hybridization has been demonstrated to contribute to the evolution of many taxa, but determining whether hybridization has anthropogenic origins, especially with human habitat modification, is a difficult problem (Allendorf et al., 2001;Dowling & Secor, 1997).Very little is currently known about the fitness of hybrids, their viability, or backcrossing direction or frequency if any occurs.Additional studies should be conducted to determine if hybridization is more widespread elsewhere around the edges of the basin where P. mcallii may overlap with P. goodei or P. platyrhinos.

| CON CLUS IONS
Molecular ecologists can quantify genetic structure across species' ranges to assess the conservation value of specific populations, but many questions remain pertaining to the practical application of results (Angert et al., 2020;Eckert et al., 2008;Hardie & Hutchings, 2010).For example, how important are rare species ecologically?How do peripheral populations persist at the range edge Smithsonian Institution (doi: 10.25572/SIHPC).We thank three anonymous reviewers for their constructive comments.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no conflict of interest.
Flat-tailed horned lizards (Phrynosoma mcallii) are distributed across the Colorado River, four states, and two nations.Their range (yellow shading) roughly corresponds to the Salton Trough.Sampling localities (n = 42 individuals) are shown as black solid circles.Red polygons denote Management Areas (MA) and Research Areas (RA) under the Rangewide Management Strategy (RMS).BB, Borrego Badlands; EM, East Mesa; OWSVRA, Ocotillo Wells State Vehicular Recreation Area; WM, West Mesa; YD, Yuma Desert; YUHA, Yuha Desert.Numbers in parentheses indicate the number of lizards included from that locality.(b) A specimen of P. mcallii from Ocotillo Wells, California.(c) Habitat near San Luis Rio Colorado, Sonora, Mexico.Photos by A.D.G.
Data overviewWe obtained a total of 81.1 million reads (an average of 1.8 million per lizard), of which 99% (80.5 million reads) passed ipyrad's initial quality filter and 92% (74.9 million reads) mapped confidently to the P. platyrhinos reference genome.Summary statistics are shown in TableS1.The average number of loci in the final assembly was 4255.Two datasets were output from ipyrad: the ingroup dataset of P. mcallii only (n = 42) and a dataset with all samples, including the outgroup and hybrids (n = 45).The ingroup SNP matrix had 20,812 variable sites with 28.9% missing sites and a sequence matrix of 942,435 bp and 27.7% missing sites.The SNP matrix containing all samples (used for phylogenetic analyses only) had 28,363 variable sites with 28.8% missing data, while the sequence matrix contained 875,955 bp with 30% missing sites.After further filtering and thinning of the ingroup dataset with VCFtools, 3767 SNPs remained for downstream population genetic and demographic analyses.
(n = 42) generated by SplitsTree shows a division(long interior 'branch')    between northwestern (green and blue labels) and southeastern (red and purple labels) populations, west and east of the Colorado River, respectively, although samples from Laguna Salada and East Mesa are intermediate (Figure2).Within the northwestern group, there is also a division between northwestern populations in the Coachella Valley and Dos Palmas (green labels) and those south and west of the Salton Sea but west of the Colorado River (blue labels).Longer terminal branches leading to the southeastern samples indicate higher genetic diversity within those populations.In the concatenated ML phylogeny (n = 45; rooted with P. platyrhinos), the two suspected hybrid individuals are deeply divergent from the remaining P. mcallii samples (Figure 3).The clade of remaining P. mcallii samples is strongly supported.Some nodes are weakly supported, with bootstrap support ≤60, and are subtended by short branches, indicating a lack of deep intraspecific genealogical lineages (Figure 3).Nodes that are strongly supported include a Coachella Valley clade, a clade comprising nearly F I G U R E 2 Unrooted phylogenetic network of concatenated RAD loci for the ingroup (n = 42) produced by SplitsTree.Colours follow Admixture results (K = 4) in Figure 4. all individuals west of the Colorado River, and a clade localized in Sonora (El Pinacate).Furthermore, northwestern individuals (especially those from Coachella Valley) are generally more deeply nested with respect to southeastern individuals in Arizona and Sonora.

Figure 4
Figure 4 shows the Admixture and DAPC results under K = 2, K = 3, and K = 4.Under K = 2, the species is separated into southeastern and northwestern populations, roughly divided by the Colorado River.At K = 3, the northwestern population is further subdivided, isolating the Coachella Valley and Dos Palmas lizards as distinct from localities south, southeast (East Mesa), and west of the Salton Sea ('southwestern' populations hereafter).Under K = 4, the El Pinacate population is separated from other populations east of the Colorado River.Although Admixture was better able to identify individuals

F
I G U R E 5 Estimated parameter surfaces under the estimated effective migration surfaces model.The left plot represents the migration parameter (m), divergence among demes, while the right plot represents the diversity parameter (q), divergence within demes.Parameter rates are colour coded on the log10 scale, with cool and warm colours representing values lower or higher (respectively) than what is expected under the null isolation-by-distance models.
The top-ranked 2D demographic model was sym_mig_size with symmetric migration and a population size change (bottleneck).The best-supported model is very strongly favoured over all others based on Delta AIC (>>>10) and Akaike Weight (1.00 vs. nearly zero for all other models).Converted demographic parameters under the sym_mig_size model and 95% credible intervals from the bootstrap analysis are presented in Table2 and Figure S3.Strikingly, these results indicate a 99% population reduction affecting both major populations (1-Nu1b/Nu1a = 0.992, 1-Nu2b/Nu2a = 0.995) that occurred very recently relative to the divergence across the Colorado River (T2/ T1 = 0.002).Under the assumed mutation rate, the 95% confidence interval for time of divergence (in years) of the northwestern and southeastern populations across the Colorado River ranges spans the late Pleistocene, 13-350 kya.The date of the bottleneck is estimated to be in the Holocene, 176-4990 years ago.

(
FigureS2) are generally consistent with the phylogenetic results in that both follow congruent hierarchies; deeper splits in the tree are reflected in the populations inferred under lower K values, and EEMS identifies corridors and barriers to migration (m) as deviations from isolation by distance.Consistent with the clustering, PCA, and phylogenetic methods mentioned earlier, our EEMS results, which make no assumptions about K values, indicate that both the Colorado River and the Salton Sea (Lake Cahuilla in prehistoric times) are barriers to migration (Figure5).Indeed, they are modelled as a single barrier to gene flow, which makes sense considering that the deepest part of the actively rifting Salton Trough captures overflow from the Colorado River, whether through natural floods and river meanderings during the Pleistocene or Holocene anthropogenic changes.
There are also interesting comparisons to be made at a local scale within the Coachella Valley.The narrow valley's topography is shaped by the San Andreas Fault, with the Little San Bernardino, San Jacinto, and Santa Rosa mountains forming steep, sharply defined escarpments on the northeast and southwest sides.This physiography creates an ecological gradient along the ~70 km length of the valley floor between the Colorado Desert in the southeast and a Mediterranean-type climate in the northwest, hinting at the potential for localized adaptation.Forty years ago, P. mcallii was widely distributed across the valley, but today occurs at only a single site (9.3 km 2 ) within the Coachella Valley National Wildlife Refuge

1.
Increase the geographic scope of the RMS to include all populations identified under K = 4, either by adding new MAs or RAs or through expanded collaboration with stakeholders in the Coachella Valley and Mexico.Angert et al. (2020), reviewing the literature on adaptation at range edges, argued that peripheral populations should not be overlooked as evolutionary dead-ends headed towards local extinction, nor should they be granted inherently superior significance solely based solely on their peripheral status.The existing design of MAs and RAs, which includes lands in California and Arizona on either side of the Colorado River, but not in Mexico, captures the deepest divergence within the species.However, the Coachella Valley and El Pinacate populations are not currently included as RAs or MAs.Fortunately, these areas are pre-existing reserves protected by the United States (USFWS) and Mexico (SEMARNAT), respectively.Although K = 2 was the best-fit model, K = 3 and K = 4 scored only marginally worse, our results and earlier discussion clearly demonstrate that these peripheral populations hold significant conservation value as repositories of genetic diversity in the existing MAs and RAs.2.Establish temporal sampling of genomic diversity for at-risk populations.We recommend that land managers establish a long-term genomic diversity monitoring programme, rangewide, to capture trends over the next century and potentially trigger direct action.As mentioned above,Vandergast et al. (2016) showed declining genetic diversity of U. inornata in the Coachella Valley between 1996 and 2008 due to anthropogenic habitat fragmentation, which occurred concurrently with extirpation of local P. mcallii populations(Barrows et al., 2021;Barrows & Allen, 2009).This was only possible through discrete temporal sampling across years, whereas our study represents individuals collected over more than two decades analysed together and establishes an important, but limited, baseline dataset for future comparisons.We recommend future studies adopt a similar (minimally invasive) temporal sampling strategy to that ofVandergast et al. (2016) for P. mcallii.As our RAD markers have been aligned to a reference genome, we suggest that researchers could use markers presented here, by either reproducing the ddRAD protocol, designing probes or primers to specific markers, or even sequencing whole genomes.If local populations become extirpated, lose diversity, show signs of inbreeding depression, or show increased population structure, managers might choose to compensate for this through genetic rescue techniques, for example, translocating juvenile lizards from one locality to another.At this time, we lack the evidence to recommend these direct interventions.3. Expand field surveys in Baja California, Mexico: More samples are needed from the Laguna Salada basin, which would require international collaboration with the U.S. RMS partners and Mexican agencies.The single sample we analysed from Laguna Salada shows admixture between the eastern and western populations, and EEMS results show this area as neutral for migration (i.e.genetic change within this basin conforms closely to the isolation-by-distance model).Suitable habitat also appears to extend further south towards San Felipe on the Gulf of California.Therefore, additional samples from Baja California have the potential to alter results and their interpretation.4. Monitor for hybrids: The known hybrids between P. mcallii and P. goodei occur in the YD MA along the bajada from the Gila Mountains west towards the Colorado River (FTHL ICC, 2003; under divergent or extreme conditions?How can managers who are contemplating translocation or genetic rescue of peripheral populations balance the opposing concerns of maximizing genetic heterogeneity through gene flow (to minimize deleterious effects of inbreeding) versus allowing for adaptation to local conditions?Our findings show that P. mcallii is an excellent candidate for pursuing these types of questions.Endemic to a desert basin formed by the rifting of tectonic plates, complete with badlands, dunes, a meandering river, and an ephemeral inland sea, P. mcallii must today cope with destroyed, degraded, or fragmented habitat, changing climate, invasive species, and increased predation.The genes of P. mcallii tell a story of late Pleistocene divergence and gene flow, northwestern range expansion, and a very severe, very recent bottleneck that yields low effective population sizes today.We cannot yet attribute the bottleneck to a single root cause, but its recency and severity point to possible anthropogenic mechanisms.This study constitutes a real-life example of how genomic data can inform applied management strategy, by identifying significant diversity in the Coachella Valley and El Pinacate reserves that fall outside of existing management and research areas, by identifying natural barriers that promote divergence, and by revealing possible mechanisms that shape intraspecific diversity.Fisher, Dan H. Foley, Lin Priest, Cameron Barrows, E. Hollenbeck, C. Itogawa, and the FTHL ICC.We thank Miguel Angel Grageda (Reserva de la Biosfera El Pinacate y Gran Desierto de Altar), R. R. Rodríguez, and the Centro Intercultural de Estudios de Desierto y Oceanos for assisting us in collecting and permits for tissue samples of both species in Mexico.California Department of Fish and Game, and Arizona Game and Fish authorized the collection of genetic samples in the United States.Laboratory work was conducted in and with the support of the Smithsonian Institution, Laboratories of Analytical Biology of the National Museum of Natural History (NMNH).Some of the analyses in this paper were run on the Smithsonian High Performance Cluster (SI/HPC), Coordinates are in WGS-84 datum.Coordinates have been rounded to a precision of two decimal points.
Colorado Desert: U. inornata in the Coachella Valley, U. notata west of the Colorado River, U. rufopunctata east of the river, and U. thurmanae of the Mohawk Dunes, Arizona, with the status of U. cowlesi of the El Pinacate region still uncertain (DeRycke