Carex specuicola is a perennial monocot endemic to the southern Colorado Plateau Desert in Arizona and Utah. The species is long-lived and occurs almost exclusively in hanging garden ecosystems (U.S. Fish and Wildlife Service, 2014), which house high numbers of endemic and rare plant species across this region (Welsh, 1989; Spence and Henderson, 1993; Fowler et al., 2007). In 1985, C. specuicola was only known from three localities (~700 individuals) and was listed as Threatened under the U.S. Endangered Species Act (ESA), with 600 square meters of critical habitat designated (U.S. Fish and Wildlife Service, 1987); since its listing, additional surveys have revealed 170 hanging garden sites occupied by C. specuicola (Rink, 2017). Occupied garden sites vary in size from 0.1–100 m² and are isolated from one another by intervening dry habitats. Vegetative spread resulting from underground rhizomes, as opposed to sexual reproduction from seed, is hypothesized to account for within-garden spread (Hermann, 1970). Rates of gene flow among gardens are unknown, though low gene flow is an assumption implicit in the species’ management (U.S. Fish and Wildlife Service, 2014) based on a lack of seed traits promoting high dispersibility, presumed pollination by wind (Friedman and Barrett, 2009), and a short duration of pollen viability (Dafni and Firmage, 2000).

Carex specuicola may be impacted by gene flow with its partially sympatric putative sister taxon, C. utahensis. Carex utahensis is a regional endemic with a more northern distribution occupying seeps, wet slopes, wet meadows, and riparian bottomlands. The species was described based on morphological differences with C. specuicola (e.g., stiff versus flexuous culm, perigynium shape, and style number; Reznicek and Murray, 2013), but morphological overlap is not uncommon, which may suggest past or current hybridization or a recent divergence time between these taxa. Evidence of gene flow between the species may have critical implications for the long-term survival of C. specuicola and its federal protection (US Fish and Wildlife Service, 2014).

We sampled leaf tissue from 190 Carex individuals across 20 sampling sites (1-12 putative independent genets per site) spanning the breadth of the C. specuicola (164 samples) distribution along with samples of C. utahensis (23 samples), C. curatorum (2 samples), and C. aurea (1 sample) ( Figure 1 and Table 1 ); the latter two species often co-occur with C. specuicola and C. utahensis in hanging garden or adjacent mesic habitat and share phylogenetic affinities, making them appropriate outgroups in phylogenomic analyses. Four sampling sites were represented by close-proximity subsites separated by 60 m to 3.27 km (1.38 ± 1.16 km), which helped us understand relatedness at small spatial scales ( Table 1 ). Proximate gardens were closely related (see Results) and were grouped for downstream analyses when indicated in the Methods. We extracted genomic DNA from leaf tissue using DNeasy 96 Plant Kits (QIAGEN, Germantown, MD, USA) following the manufacturer’s protocol. We then prepared ddRADseq libraries with individually barcoded samples using EcoRI and MspI restriction enzymes (Peterson et al., 2012). We pooled individual libraries and amplified them via polymerase chain reactions. We then isolated amplified genomic fragments ranging in length from 400 to 600 base pairs (bps) with a Pippin Prep (Sage Science, Beverly, MA, USA). We sequenced the final pooled library on a HiSeq 4000 (Illumina, San Diego, CA, USA) at the University of Oregon Genomics and Cell Characterization Core Facility to generate single-end, 100 bp reads.

After removing five samples that yielded low quality DNA ( Table 1 ), we processed and aligned raw sequence data of 185 samples using Stacks 2.41 (Rochette et al., 2019). First, we cleaned raw data by removing reads with uncalled bases or an average phred-scaled quality score below 22 within a 15 bp sliding window. We also trimmed the last 5 bp because quality scores decline sharply towards the end of reads. We determined the appropriate de novo assembly and filtering parameters in Stacks following Rochette and Catchen (2017). We first processed a test dataset of 19 samples (16 C. specuicola, 3 C. utahensis) while varying the parameters M (the maximum nucleotide distance allowed between stacks) and n (the number of mismatches allowed between loci during catalog construction) from 1 to 9 and while keeping M = n and setting m = 3 (minimum depth of coverage to create a stack). Following this procedure, we processed the full dataset with M = 4 and n = 4 because the total number of loci retained in 80% of individuals (r80) and the distribution of the number of single nucleotide polymorphisms (SNPs) per locus appeared to stabilize at these values ( Figure S1 ). We calculated genetic summary statistics across sampling sites of C. specuicola and C. utahensis in the populations program in Stacks using filtering parameters that removed loci with heterozygosity > 75% (–max-obs-et 0.75) and loci not present in at least 50% of individuals in a sampling site (-r 0.5).

To ensure that downstream analyses were not influenced by admixture with C. utahensis, we inferred relationships among Carex species using a maximum-likelihood phylogenetic analysis in PhyML 3.0 (Guindon et al., 2010). Model selection was implemented using smart model selection (Lefort et al., 2017) and nodal support values were calculated from 100 bootstrap iterations. The phylogeny was rooted using C. aurea and C. curatorum as outgroups. Based on the placement of Chinle Wash between C. specuicola and C. utahensis (see Results), we performed f₃ -admixture tests with treemix v1.13 (Pickrell and Pritchard, 2012), which use information on the frequency of shared SNPs among three populations to infer whether the relationship between the three populations is best explained by admixture as opposed to a simple bifurcating tree (Peter, 2016). A negative f₃ statistic indicates that a population results from admixture. We calculated f₃ statistics for the Chinle Wash (C. specuicola) site relative to parental C. specuicola and C. utahensis sites using a black-jackknife approach with 160 bins of 100 random SNPs to estimate standard errors.

To gain a better understanding of the processes influencing hanging gardens, we estimated mean pairwise Hudson’s F_ST with KRIS (sensu Bhatia et al., 2013) between hanging gardens with hierfstat (Goudet and Jombart, 2020) in R 4.1.2 (R Core Team, 2021). These data were leveraged to test for isolation-by-distance (IBD) patterns across the species’ range and for each genetically defined group identified by Structure analyses (see below). In all cases, we used Mantel tests (1x10⁵ permutations) and linear regressions to test for a relationship between genetic and geographic distance matrices (Mantel, 1967; Bougeard and Dray, 2018). Next, we calculated mean inbreeding coefficients (F_IS ) in hierfstat and mean relatedness (unadjusted mean A_jk ; Yang et al., 2010) between samples within sites using the -relatedness function of vcftools (Danecek et al., 2011).

To investigate intraspecific genomic structure, we employed two approaches including 1) Bayesian clustering implemented in Structure v2.3.4 (Pritchard et al., 2000) and 2) principal coordinates analysis (PCA) implemented in the adgenet package (Jombart, 2008) in R. Given that these methods complement one another due to their reliance on different assumptions (Jombart et al., 2008; Frantz et al., 2009), we assess their results in tandem to determine biologically relevant, genetically defined groups instead of using model selection methods that can be misleading (i.e., Janes et al., 2017). For Structure, we tested K values (representing the number of genetically defined populations) from 1 to 8, performing 5 iterations under each K value with 20,000 burn-in and 100,000 regular repetitions. We excluded Chinle Wash individuals because they showed evidence of C. utahensis genomic identity (see Results).

To further investigate the spatial distribution of genetic diversity, we used the program EEMS (Estimated Effective Migration Surfaces; Petkova et al., 2016), which models the relationship between genetic variation and geography using Markov chain Monte Carlo sampling to explore large parameter spaces (Petkova et al., 2016). This procedure results in a visualization illustrating the geographic areas where occupied sites are expected to have higher or lower effective diversities, on average, based on patterns across the entire empirical data. We ran EEMS assuming a stepping-stone migration model with either 50, 100, or 200 subpopulations (termed demes in EEMS) spaced in a grid across a triangular range encompassing all sampled localities. We excluded Moenkopi Wash because we sampled only one individual from this site and grouped proximate gardens (see Table 1 ). We performed 3 independent runs for each subpopulation value to ensure model convergence. We used default hyper-parameter values and tuned the proposal variances until proposal acceptance rates reached 20-40%. We ran EEMS for 2 million iterations with a burn-in of 1 million iterations and thinning iteration of 9999.

We performed demographic modeling to estimate divergence times (t), ancestral effective population sizes (N_anc ), contemporary effective population sizes (N_e ), and levels of historical gene flow between pairs of closely related C. specuicola hanging garden sites. Specifically, we explored the pairwise evolutionary history of the Forty Caves and Tse Ya Toe Spring sites; Geshi and Canyon del Muerto sites; and Cigarette Spring and Road Canyon sites. Each pair of hanging gardens was closely related according to PCA, Structure, and phylogenetic analyses, though we excluded two Geshi individuals that may have resulted from a recent dispersal event (see Results).

We used ∂a∂i (Gutenkunst et al., 2009) to infer demographic parameters using a folded site frequency spectrum (SFS; i.e., major/minor allele calls rather than derived/ancestral). We scaled parameters by θ, assuming the tomato neutral mutation rate (m = 1 × 10^-8; Lin et al., 2014) and a 337:1 ratio of callable sites:SNPs, based on output from Stacks. For each pair of sites, we tested two demographic models: 1) no migration, 2) isolation-with-migration (IM). The ‘no migration’ model involved a split between sites at t generations in the past and allowed for instantaneous population size changes immediately following divergence. The ‘IM’ model was identical but included a symmetrical migration rate (m; migrants per generation) between sites immediately following divergence t generations ago that remained constant up to the present day.

For each model, we performed 20 independent runs starting with parameter values sampled randomly across the following uniform distributions: 0.01 < n < 100.00, 0 < 2N_anct < 20, 0 < 2N_ancm < 5; where n is the relative change in population size compared to N_anc . We determined which of the two maximum likelihood demographic models produced the best overall fit to the data using a composite-likelihood ratio test with the Godambe Information Matrix (GIM; Coffman et al., 2016) and Akaike Information Criterion scores (AIC; Akaike, 1973). We estimated confidence intervals for each parameter with the Godambe Information Matrix with 100 bootstrap data sets comprising 50% of SNPs randomly selected from the full data set. We validated the final model by comparing the predicted SFS to the observed SFS for each population.

The maximum likelihood phylogeny of the four Carex species demonstrates the close and highly supported relationships within and between C. specuicola and C. utahensis ( Figure 2 ). Most nodes have support values >95%, except for a few instances within the C. specuicola and C. utahensis clades. All proximate garden sites (i.e., those located within 3.27 km) clustered together with high bootstrap support. The only exception to the pattern of close intraspecific relationships is represented by Chinle Wash, which is sister to a clade containing the remaining C. specuicola gardens. As this may be a sign of past or ongoing gene flow with C. utahensis, we calculated f₃ statistics to identify a signal of hybridization. However, we found no C. specuicola × C. utahensis comparison that yielded a negative f₃ value ( x ¯ ± SEM = 0.0202 ± 0.0018).

Observed and expected heterozygosities within C. specuicola were generally low ( x ¯ = 0.04 and 0.03, respectively), as were nucleotide diversities ( x ¯ = 0.04; Table 1 ). Private alleles among sites ranged from 27-418 ( Table 1 ). Genetic differentiation (i.e., F_ST ) between pairs of hanging garden sites was high, with an overall x ¯ ± SD = 0.49 ± 0.14 ( Table S1 ), indicating substantial population structure even across small spatial scales. For instance, Adahchiiyahi Wash and A-Canyon are separated by only ca. 4 km but had a F_ST = 0.4. Across 325 pairwise F_ST values, the highest F_ST was between Tse Ya Toe Spring and Tes Nez Iah (F_ST = 0.68). Pairwise F_ST values between subsites ranged from -0.10–0.38 ( Table S1 ), with both the highest and lowest F_ST values among Walker Creek subsites ( Table S1 ). Geographic and genetic distance were significantly correlated in both Mantel tests and linear regressions across all hanging gardens and within Structure defined genetic clusters ( Figure 3 ; Table S2 ).

We estimated levels of inbreeding (F_IS ) and relatedness (A_jk ) to investigate evidence of clonal reproduction within hanging gardens. F_IS and A_jk were negative at most sites, suggesting that neither inbreeding nor clonal reproduction is widespread ( Table 1 ). However, we found significantly positive F_IS values across two sites (Forty Caves and Geshi), which may indicate recent colonization, inbreeding, or cryptic population structure. At Geshi, F ¯ I S was no longer strongly elevated when calculated separately for individuals that group with Inscription House Spring ( F ¯ I S =   −0.14) or Canyon del Muerto ( F ¯ I S = −0.002), supporting cryptic population structure as the probable driver of the pattern.

Structure and PCA analyses resolve similar patterns of genetic structure across C. specuicola hanging gardens. At K = 2, hanging gardens are composed entirely of one genetic axis or the other ( Figure S2 ) and only two sampling sites (Forty Caves and Moenkopi Wash) show evidence of admixture. At K = 3, one of the genetic axes splits into two, while the other genetic axis remains relatively consistent. The K = 3 Structure solution reiterates relationships discerned by the PCA, in which individuals group into three primary clusters ( Figure 1B ). All individuals from a garden cluster with one another in PC space except for four, including two from Geshi that are more closely related to Inscription House/Binne Etteni, one from Jackass that clusters with Shonto, and one from Cigarette Spring that approaches individuals from Rock Point ( Figure 1C ). In relation to geography, genetic ancestry of individuals from the Sheik and Canyon del Muerto sites was more closely aligned to individuals from distant sites (e.g., Adahchiiyahi Wash and Geshi, respectively) than to individuals from the nearest sites. We also observed mixed ancestry at the Snow Flat Spring Cave, Rock Point, and Moenkopi Wash in which the minor ancestry component was present in geographically distant gardens ( Figure 1B ). The K = 4 model highlights a close relationship between Geshi and Canyon del Muerto that is not evident in the PCA, as well as illustrates the mismatch of the two aforementioned Geshi individuals ( Figure S2 ). Structure analyses at higher K values did not provide biologically relevant or interpretable patterns.

EEMS runs produced strong correlations between observed and expected genetic dissimilarity both within and between demes, indicating good model fit (Petkova et al., 2016). Effective diversity surfaces were similar across different deme models ( Figures 4 , S3 ); as a result, we report the surface developed using 200 subpopulations because it most closely reflects the number of known occupied hanging garden sites (170 sites; Rink, 2017). The effective diversity surfaces revealed relatively high genetic diversity (25-30% higher than average) in the far western portion of the range (i.e., Binne Etteni, Inscription House Spring; Geshi; Forty Caves; and Tse Ya Toe Spring; Figure 1 ) and lower than average diversity to the east. This is generally in agreement with the expected heterozygosity and nucleotide diversity of sites, which were both higher in the western portion of the range ( Table 1 ).

Modeling of population history between pairs of closely related hanging garden sites generally revealed small effective population sizes and recent divergence times, although support for gene flow between sites was mixed. We found strong support for an IM model over a no migration model (GIM p< 0.001, ΔAIC = 40) between Forty Caves and Tse Ya Toe Spring (F_ST = 0.4, 4 km distance); however, confidence intervals for divergence time, migration rate, and population size were wide ( Table 2 ). Our maximum likelihood estimates place the Forty Caves–Tse Ya Toe Spring divergence time at 2595 generations ago with ca. 2.5 migrants per generation following divergence. Forty Caves had the largest effective population size of any hanging garden sites (maximum likelihood N_e = 1089), while Tse Ya Toe Spring was characterized by a very small effective population size (maximum likelihood N_e = 48). The common ancestor of Forty Caves and Tse Ya Toe Spring also had a relatively large effective population size of an estimated 4340 but confidence intervals around these parameter estimates were large. Between Canyon del Muerto and Geshi (F_ST =0.26, 127 km), we found support for a no migration (ΔAIC = 7) model that included very low effective population sizes (Canyon del Muerto N_e = 8–17; Geshi N_e = 0–764) and recent divergence between sites (334–714 generations ago; Table 2 ). Finally, for Cigarette Spring and Road Canyon (F_ST = 0.03, 3 km distance), no migration and IM models explained variation in the data equally well (ΔAIC = 1.7). Effective population sizes were small (e.g., Cigarette Spring N_e = 12-261; Road Canyon N_e = 0-164) and divergence times were recent (0-302 generations ago; Table 2 ) regardless of whether gene flow was included in the model.

Carex specuicola and C. utahensis are closely related Colorado Plateau endemics within the “Carex parryana” complex of Carex sect. Racemosae (Reznicek and Murray, 2013; Massatti et al., 2016). Our analyses suggest that C. specuicola and C. utahensis form reciprocally monophyletic groups ( Figure 2 ) and that hybridization appears to be rare. The conflicting results between the f₃ tests and the phylogenetic hypothesis for Chinle Wash are interesting, though a non-negative f₃ statistic does not mean that admixture did not occur; strong genetic drift, which appears common across C. specuicola hanging gardens, can quickly erode signatures of admixture in the f₃ statistic (Patterson et al., 2012). Thus, the lack of an f₃ admixture signature in Chinle Wash may be consistent with relatively old admixture. Investigating the broader genomic and evolutionary consequences of hybridization and dating the divergence between these species are intriguing areas of future research; such studies may also provide additional insight into the historical factors influencing diversification across this spatially complex and environmentally dynamic region (Kadereit and Abbott, 2021).

We explored patterns of population structure and demographic history to infer the factors influencing the formation of and interaction among C. specuicola hanging gardens. We uncovered a broad spatial pattern of genetic differentiation consistent with IBD ( Figure 3 ) and patterns of ancestry wherein geographically adjacent sites tended to be closely related ( Figure 1 ). At small spatial scales, as represented by our sampling of proximate subsites, individuals were almost always found to be highly related ( Table S1 ). Moreover, demographic inference supported a model of continued gene flow since divergence between Forty Caves and Tse Ya Toe Springs (separated by 3.8 km; Table 2 ), and we found that four individuals did not cluster with individuals from the same garden, but rather with individuals at nearby gardens, potentially indicating recent dispersal events ( Figure 1C ). Analyses supporting gene flow across the species’ range may explain why inbreeding coefficients and relatedness values are low despite F_ST values that may indicate strong genetic drift.

In tandem with evidence indicating that gene flow among adjacent sites influences patterns of genetic structure, we uncovered genetic patterns indicating C. specuicola is capable of long-distance dispersal events. For example, although the Geshi and Canyon del Muerto sites are 127 km apart, they were found to be closely related based on PCA, Structure, and phylogenetic results. Moreover, their estimated divergence time (519 generations) is more recent than is the divergence time of the spatially proximate Forty Caves and Tse Ya Toe Springs sites (2595 generations; Table 2 ). Furthermore, we found that Sheik individuals are more closely related to individuals from Adahchiiyahi Wash and A-Canyon (72 and 76 km, respectively) than to nearby sites like Road Canyon (21 km; Figure 2 ). Long-distance dispersal requires a mechanism, which is not evident in the intrinsic seed traits of Carex species; for example, seeds lack structures that facilitate dispersal by wind (e.g., a pappus as in Asteraceae species) or by hitchhiking (e.g., barbed seeds as in some Boraginaceae species). However, Wilson et al. (2008) found that Carex seeds can pass intact and viable through the digestive systems of waterfowl and gallinaceous birds, potentially implicating zoochory. Within desert ecosystems, mesic plant communities like hanging gardens may be important sources of food, water, and shelter for birds (Malanson and Kay, 1980; Welsh, 1989), thereby providing a direct avenue of propagule dispersal among the gardens distributed across otherwise heterogeneous, complex, environmentally unsuitable terrain.

Population genetic patterns and inferred demographic processes suggest that C. specuicola may have metapopulation-like dynamics. Metapopulations are governed by recurrent extinction and colonization processes (Pannell and Charlesworth, 2000), leading to temporary occupancy of habitat patches (Thomas, 1994). In this system, the recent estimated divergence times between pairs of hanging garden sites ( Table 2 ) are consistent with rapid population turnover of occupied/unoccupied sites; furthermore, the range of C. specuicola’s hanging garden F_ST values reflects biological systems in which colonization occurs regularly (Knowles and Alvarado-Serrano, 2010). Population census data support these processes; for instance, Hevron (1990) recorded negative survey results (i.e., no C. specuicola found) at 10 sites in the Forty Caves area that were subsequently resurveyed by Rink (2017), who located C. specuicola at 6 of the 10 sites. Likewise, extirpation at one site has been documented (Rink and Hazelton, 2014). Extinction and colonization within metapopulations is expected to strongly reduce population genetic diversity (i.e., effective population size; Harrison and Hastings, 1996; Pannell and Charlesworth, 2000). Consistent with this, we find low relative genetic diversity and small effective populations across most of the range ( Figure 4 , Table 1 ). However, the western portion of C. specuicola’s range has relatively high genetic diversity ( Figure 4 ), which may reflect larger and more stable populations and thus a source of the metapopulation; this hypothesis requires further testing with additional census and genetic data. In light of geographic patterns of genetic diversity, hanging garden ecosystems may be biogeographically similar to oceanic and sky island systems, where species composition is determined by dispersal limitations, extirpation, and habitat size (MacArthur and Wilson, 1963; Malanson and Kay, 1980; Atwood et al., 2011; Vásquez et al., 2016; He et al., 2019).

We presume that similar intraspecific dynamics occurred throughout the phylogeographic history of C. specuicola, specifically regarding gene flow and dispersal. Considering the paleorefugium hypothesis (our null hypothesis), our results do not support that this system is dominated by persistent extirpations (Nekola, 1999) to the point that the species’ long-term survival is at risk (assuming the maintenance of suitable hanging garden habitat). In a scenario where C. specuicola was widespread across suitable habitat during cooler and wetter periods (Spence and Henderson, 1993; Coats et al., 2008) and maintained similar levels of dispersal, population genetic structure was presumably more homogenous across the species’ distribution. Subsequent environmental warming and drying restricting C. specuicola to paleorefugia (i.e., hanging gardens) would have restricted gene flow, but it is unclear whether range fragmentation in the face of ongoing gene flow would be able to produce the patterns of small effective population sizes and large F_ST values inferred herein.

The alternative hypothesis to a system dominated by paleorefugia – that C. specuicola has always been a hanging garden specialist – seems more consistent with our data. Namely, groups of individuals within gardens do not display large inbreeding coefficients or a lack of dispersal. While support for the paleorefugium hypothesis in a co-distributed hanging garden species (Anticlea elegans subsp. vaginata – a monocot in the family Melanthiaceae) was in part garnered from the distribution of the species’ sister taxon in montane habitat of the San Juan Mountains to the northeast of this study region (Palmquist et al., 2015), we note that C. specuicola is within a clade of species that is not restricted to higher elevations (Massatti et al., 2016). In other words, while vicariance presumably influenced diversification of Anticlea elegans as the climate warmed after the Last Glacial Maximum and the parent taxon retreated to (or only persisted in) higher elevations, a similar diversification event involving a montane sister taxon is not supported for C. specuicola. However, this hypothesis does not preclude C. utahensis (or an ancestor of the two species) having a more widespread distribution during cooler and wetter times, and perhaps a drying and warming climate may have reduced gene flow between these species to the point where they are now largely reproductively isolated. Similar diversification processes are likely responsible for diversification within Carex (Martín-Bravo et al., 2019, Hodel et al., 2021), as well as within other clades distributed across topographically and environmentally heterogeneous landscapes (Jones et al., 2021; Kadereit and Abbott, 2021). A spatially and temporally explicit modeling framework may be the best way to test alternative hypotheses (e.g., Bemmels et al., 2016; Massatti and Knowles, 2016; Ortego and Knowles, 2020). Regardless of C. specuicola’s true history, our data indicate that the species maintains gene flow and can colonize hanging gardens through dispersal processes.

Understanding dispersal, colonization, gene flow, and interspecific dynamics in rare species distributed across isolated habitat is essential for developing effective conservation strategies. Our data indicate that hybridization with C. utahensis is rare and geographically restricted, which alleviates admixture concerns regarding the persistence of C. specuicola. Furthermore, gene flow and dispersal processes influence population structure and genetic diversity. In light of these results, we suggest that future habitat loss due to climate change, drought, or other factors represent the major threat to this species, given that fewer hanging gardens translates to less habitat for C. specuicola to potentially inhabit. As hanging garden habitat is hard to survey given accessibility issues, a management plan involving recurrent genetic surveys may indicate whether declining habitat quality is influencing genetic diversity and structure (Paz-Vinas et al., 2018). To this end, we note that our investigation provides baseline data that can inform future surveys. The preservation of multiple, genetically distinct populations is crucial to maintain existing genetic diversity and to support colonization of new sites if hanging garden populations stochastically extirpate. Further, the maintenance of redundant sites with similar genetic profiles may prevent stochastic events from reducing genetic richness. Given the congruence of genetic patterns between PCA and the Structure K = 3 analysis, we identify three primary populations that may have utility in management applications, including the East population (represented herein by Road Canyon, Cigarette Springs, Walker Creek sites, Rock Point, and Tes Nez Iah); the Central population (represented herein by Sheik, Adahchiiyahi Wash, A-Canyon, Jackass, Peach Spring, and Shonto); and the West population (represented herein by Geshi, Binne Etteni, Inscription House, Moenkopi Wash, and Canyon del Muerto). However, due to the dispersal biology of the species, we emphasize that it may be difficult to predict the genetic ancestry of a site with high confidence (e.g., Massatti and Winkler, 2022), especially as distance increases from a site with known ancestry.

Although our work provides important insights into the biology of this species and potential conservation strategies, future research is essential to better understand whether it can persist under current conditions and how it may respond to future challenges. In particular, characterizing dispersal distance and the spatial structure of suitable habitat are crucial for understanding capacity for long-term persistence of metapopulations (Grilli et al., 2015; Hanski et al., 2017). In addition, investigating potential adaptive loci could provide insights into genetic or environmental factors that influence metapopulation size and persistence (e.g., PGI in butterflies, Hanski et al., 2017). For instance, future work could explore spatial patterns and potential genetic underpinnings of drought tolerance, which could be integrated into conservation strategies to facilitate adaptation under climate change. Ultimately, strategies for managing rare species depend upon understanding how historical and contemporary processes affect intraspecific patterns and processes with an eye on what the future may hold; genetic investigations are a prime mechanism to meet these needs.