Genomics‐informed delineation of conservation units in a desert amphibian

Abstract Delineating conservation units (CUs, e.g., evolutionarily significant units, ESUs, and management units, MUs) is critical to the recovery of declining species because CUs inform both listing status and management actions. Genomic data have strengths and limitations in informing CU delineation and related management questions in natural systems. We illustrate the value of using genomic data in combination with landscape, dispersal and occupancy data to inform CU delineation in Nevada populations of the Great Basin Distinct Population Segment of the Columbia spotted frog (Rana luteiventris). R. luteiventris occupies naturally fragmented aquatic habitats in this xeric region, but beaver removal, climate change and other factors have put many of these populations at high risk of extirpation without management intervention. We addressed three objectives: (i) assessing support for ESUs within Nevada; (ii) evaluating and revising, if warranted, the current delineation of MUs; and (iii) evaluating genetic diversity, effective population size, adaptive differentiation and functional connectivity to inform ongoing management actions. We found little support for ESUs within Nevada but did identify potential revisions to MUs based on unique landscape drivers of connectivity that distinguish these desert populations from those in the northern portion of the species range. Effective sizes were uniformly small, with low genetic diversity and weak signatures of adaptive differentiation. Our findings suggest that management actions, including translocations and genetic rescue, might be warranted. Our study illustrates how a carefully planned genetic study, designed to address priority management goals that include CU delineation, can provide multiple insights to inform conservation action.


| INTRODUC TI ON
Appropriate identification of conservation units (CUs) is essential for effective management because CUs inform conservation listing status (e.g., under state, provincial, or federal jurisdictions) and management options available to support population recovery and persistence. Conservation unit is a general term for an array of biological and legal categories used to define subspecific population units for the protection and management of intraspecific diversity.
Because inappropriate delineation of CUs can waste conservation dollars, hamper recovery efforts and result in inappropriate management actions (Weeks et al., 2016), careful evaluation of CU delineation is crucial for species of conservation concern.
One of the most commonly used CUs is the evolutionarily significant unit (ESU), which represents major components of intraspecific variability based on genetic, phenotypic, ecological, geographical and life history characteristics. There are many definitions and sets of criteria for defining ESUs (reviewed in Fraser & Bernatchez, 2001), but here we will focus on the parameters established by Waples (1991): ESUs show substantial reproductive isolation from other conspecific populations and represent an important component of the evolutionary legacy of the species. We use this definition because it best reflects the analogous policy unit under the U.S. Endangered Species Act (ESA), the "distinct population segment" (DPS). Note that DPS designation is a policy decision made by the agencies that implement the ESA, so this study will not evaluate DPSs. Management units (MUs), which are often nested within ESUs, are demographically independent, meaning that population dynamics rely on local birth and death rates rather than immigration (Palsbøll et al., 2007). MUs typically warrant separate management due to their importance for maintaining functional ecological and evolutionary processes within the larger ESU (Moritz, 1999). The conservation of multiple MUs within an ESU is considered essential for maintaining long-term persistence, especially in species showing metapopulation dynamics (Gustafson et al., 2007;Weckworth et al., 2018).
These biologically defined CUs interface with policy definitions in diverse ways across jurisdictions. For example, under the ESA, legal protections can be extended not only to species and subspecies, but also DPSs of vertebrates. Since the ESA does not specify the basis for DPSs, the two agencies responsible for implementing the ESA have defined the term based on two main axes: "discreteness" (reproductive isolation) and "significance" (substantial contribution to evolutionary legacy and potential for persistence), similar to the Waples (1991 ESU definition above (USFWS & NMFS, 1996;Waples, 2006). Despite this clarification, CU delineation is often a complex task under the ESA, with ramifications for candidate listing determinations and extinction risk assessments (Waples & Lindley, 2018).
Genetic data can play an important role in informing the complexities of CU delineation and subsequent management and policy decision-making. Neutral genetic markers, including mitochondrial DNA, microsatellites and single nucleotide polymorphisms (SNPs), can be used to evaluate reproductive isolation in ESUs and, to a limited extent, demographic independence in MUs Lowe & Allendorf, 2010;Palsbøll et al., 2007).
Additionally, large genome-wide data sets have improved our ability to identify potential adaptive differentiation, informing the "evolutionary legacy" criterion in ESUs. Importantly, however, it is not always clear how to delineate CUs based solely on genetic or genomic data. For example, there is no single, simple relationship between demographic independence and genetically informed estimates of population structure (Palsbøll et al., 2007), which means that genetic data alone can often provide little information on the demographic independence criterion for MUs (Lowe & Allendorf, 2010). These limitations are important to consider in the context of management objectives for species, given the importance of CUs for listing and conservation decision-making (Taylor & Dizon, 1999).
We provide a case study in the use of genomic data to inform CU delineations in a desert amphibian. By combining multiple genomic analyses with environmental data and information on dispersal capacity and site occupancy, we evaluate evidence for ESUs and MUs in Nevada populations of the Columbia spotted frog (Rana luteiventris). R. luteiventris occupies a large range that has been divided into three major clades (Funk et al., 2008; Figure S1). We analysed genome-wide SNP data from sites throughout Nevada, representing most of the range of the Great Basin clade of R. luteiventris and all three occupied regions in Nevada: the Jarbidge-Independence (the core part of the range), the nearby but isolated Ruby Mountains and the geographically disjunct Toiyabe (Figure 1). While these populations occupy naturally fragmented aquatic habitats in this xeric region, anthropogenic activities over the past 200 years, such as beaver extirpation, stream incision, water diversions, wetland draining, agricultural irrigation and urbanization, have substantially altered the regional hydrology. These activities, along with changes in climate, have probably eliminated or altered many of the cool lentic and lotic habitats preferred by R. luteiventris Welch & MacMahon, 2005). Across the Great Basin, extant populations of R. luteiventris persist in locations where temperature and precipitation have changed the least in the last 100 years, suggesting that remaining populations may be more isolated than they were historically . Further evidence for this isolation effect was found in a survey of 220 sites across the Great Basin, which reported that watershed occupancy increased with precipitation and decreased with temperature (Smith & Goldberg, 2020).
The compounding impacts of invasive species (i.e., American bullfrog [Lithobates catesbeianus], predatory fishes), disease, drought and climate change places many of these populations at high risk of extirpation without management intervention Mims et al., 2020;Pilliod et al., , 2021.
The U.S. Fish and Wildlife Service (USFWS) first identified the Great Basin clade of R. luteiventris as a DPS eligible for listing under the ESA in 1997. At that time, the DPS was placed on the candidate list, meaning there was sufficient information to propose the DPS as threatened, but the proposal was precluded by higher priority listing activities. The Great Basin DPS remained on the candidate list until a "not warranted" ruling resulted in its removal in 2015 (USFWS, 1993(USFWS, , 2015a. Candidate listing removal was due in large part to proactive conservation actions by diverse stakeholders in the form of ongoing 10-year Conservation Agreement and Strategy plans (Columbia Spotted Frog Technical Team, 2003a, 2003b, 2015McAdoo & Mellison, 2017). Despite these efforts, however, many R. luteiventris populations in the Great Basin continue to experience declines (USFWS, 2015b Figure 1). Additionally, genetic data are expected to inform other priorities, such as identifying dispersal corridors, evaluating genetic diversity within populations and adaptive differentiation among populations, and identifying suitable populations for actions such as genetic rescue (Bell et al., 2019;Frankham, 2015;Whiteley et al., 2015). These actions are especially important for the longterm viability of patchily distributed species such as Great Basin R. luteiventris, which have relatively low dispersal capacity (~5-6 km observed: Funk, Greene, et al., 2005;Reaser, 1996) and small population sizes, putting them at risk of increased genetic drift, inbreeding, decreased fitness and extirpation.
To inform ongoing management of R. luteiventris in Nevada, we addressed three primary research questions: (i) Are there ESUs F I G U R E 1 Locations of 31 Rana luteiventris sites sampled for genomic data in Nevada (large shapes with four-letter bold labels). Small grey points show other known occupied sites. Sampled sites are colour coded based on watershed/current management unit (labelled in italics). Major geographical regions are labelled in bold. The Toiyabe region (inset) is disjunct from the remainder of the range and represents the southernmost locations in the species range. Extent map shows the Jarbidge-Independence and Ruby Mountains regions in red, with the Toiyabe region in black; grey range map from the USGS Gap Analysis Project (2018).
within Nevada populations of the currently defined Great Basin DPS? (ii) Do genetic data and connectivity analyses support the current delineation of MUs in Nevada? If not, how should MUs be revised to better support management goals? (iii) How should MUs, sites and movement corridors be managed and prioritized for conservation actions such as genetic rescue and reintroduction efforts? 2 | MATERIAL S AND ME THODS

| Sample collection
Toe and tail clips were collected between 2010 and 2019 from 31 sites ( Figure 1). Individuals were captured at breeding ponds and released at their capture location. Sampling was conducted under authorization of USFWS (Reno, NV) and Nevada Department of Wildlife (Elko, NV).

| Genotyping
We used Qiagen DNeasy Blood and Tissue Kits for DNA extraction, adding 4 μl of RNase after tissue digestion. We quantified DNA using Qubit dsDNA assays and verified quality in a subset of samples using agarose gel electrophoresis. We used double digest restriction-site associated sequencing (ddRADseq) to develop genome-wide SNPs (Peterson et al., 2012), testing five combinations of restriction enzymes to identify the combination producing fragment sizes in the target range (250-450 bp). Using NspI and EcoRI we digested ~1 μg of DNA per individual. For library production, we used a Pippin Prep size selection of 346-456 bp, KAPA HiFi Hotstart Ready Mix for polymerase chain reaction (PCR) amplification, degenerate barcodes on the P2 adapter to allow for removal of PCR duplicates (Schweyen et al., 2014), and Dynabeads cleanup to retain molecules with P1 and P2 adapters. We sequenced one library with the University of Oregon Genomics Core Facility (HiSeq-4000, 100 bp pairedend reads), and the remaining nine libraries with Novogene Corp.
(HiSeq-4000, 150 bp paired-end reads). The first two libraries contained 24 and 36 individuals; all remaining libraries contained 40 individuals. Individuals from sampled locations were split across libraries to avoid confounding location and library effects.
We used fastqc v.0.11.8 (Andrews, 2019) to assess data quality. We removed PCR duplicates using clone_filter from stacks version 2.41 (Rochette et al., 2019). Using trim galore! version 0.6.4 (Krueger, 2019) and cutadapt version 2.5 (Martin, 2011), we removed degenerative barcodes and applied quality screening using a Phred cutoff = 20 and auto-detect adapter screening (stringency = 6). We used stacks process_radtags to demultiplex, trim reads to 90 bp, remove reads with an uncalled base or low-quality scores, and rescue barcodes and cutsites with at most one mismatch.
We used within-and across-library replicates to optimize three de novo stacks parameters: minimum stack depth (-m), distance between stacks (-M) and distance between catalog loci (-n). We used 19 individuals (each with a replicate sample) from 16 sites, covering the geographical range of sampling, and tested 13 combinations of parameters starting from default settings (-m 3, -M 2, -n 1). We evaluated mean coverage, number of assembled loci, number of polymorphic loci shared across 80% of individuals, number of SNPs (Paris et al., 2017) and population/duplicate clustering using multidimensional scaling plots (Gugger et al., 2018) in plink version 1.90 (Chang et al., 2015).
We applied a coarse filter using populations: -p 1, -R 0.3, --min_mac 2. We discarded four individuals with high missingness (>55%) before filtering with the r package radiator (Gosselin et al., 2020) using the following settings: retain common markers across 31 sites; retain SNPs with global minor allele count (MAC) >9; retain loci with coverage between 6 and 100, genotyped in at least 80% of individuals, and with ≤7 SNPs per locus; retain one SNP per locus, keeping the SNP with highest MAC. We removed 204 SNPs that deviated significantly (p < .05) from Hardy-Weinberg proportions in five or more sites (using sites with ≥10 individuals). We removed 1328 SNPs that were heterozygous in all individuals (potential paralogues) in at least one site (using sites with seven or more individuals genotyped). We removed 476 SNPs associated with a weak lane effect in the first library by identifying the contributing SNPs using redundancy analysis (RDA) from the package vegan version 2.5-6 (Oksanen et al., 2019). We first imputed missing values using snmf in lea version 3.1.2 (Frichot & François, 2015), testing K = 1-15 and alpha (regularization parameter) = 10, 100 and 1000, with 25 repetitions, 200 iterations and 5% cross-entropy withholding (final settings: K = 10, alpha = 10). The RDA used the imputed data set as the response and a dummy variable for library as the predictor. Using a threshold of ±4 SD from the mean loading on the constrained axis, we identified the loci contributing to the lane effect and removed them ( Figure S2).
We used the whoa package version 0.0.2.999 (Anderson, 2019) to assess heterozygote miscall rates for all sites with sample sizes ≥11, using polymorphic SNPs genotyped in at least 67%-86% of individuals (depending on sample size). We did not remove putative siblings based on recommendations from Waples and Anderson (2017).
We used r versions 3.6.3 and 4.1.0 (gravity models only, R Core Team, 2020).

| Objective 1: ESUs
We identified candidate adaptive markers using partial RDA (pRDA) at the site level (Forester et al., 2018), a genotype-environment association analysis that identifies genomic signatures of selection related to environmental variables. We removed low-frequency SNPs (variable in <3/31 sites, retaining 38,712 SNPs) and used three uncorrelated (|r| < 0.6), biologically relevant environmental predictors derived from the 800-m PRISM climate normal (1981( , Daly et al., 2008: average summer precipitation (June-August, when frogs are at breeding sites or moving among foraging areas), average autumn precipitation (September-November, when frogs are moving to overwintering sites) and average winter minimum temperatures (December-February, a measure of winter severity). We accounted for population structure using PC axes (Capblancq & Forester, 2021). Although the principal components analysis (PCA) screeplot ( Figure S3) showed a break in variance at PC1, we retained seven PCs based on population structure results (below). We used a candidate threshold of loci loading ±3.5 SD from the mean loading of the first two constrained axes. We used a site-based RDA of candidate markers to identify drivers of adaptive divergence and an individual-based PCA to investigate adaptive differentiation.
We assessed reproductive isolation by evaluating population structure using neutral markers (genetic data with candidate adaptive SNPs removed). We calculated pairwise F ST (Weir & Cockerham, 1984) using stampp version 1.6.3 (Pembleton et al., 2013) and visualized population structure of the individual-based data using a model-free approach, PCA, and an ancestry estimation method, admixture version 1.3 (Alexander et al., 2009). For PCA, we used the snmf-imputed data. For admixture, we tested K = 2-20 and used a cross-validation plot to determine the best K(s).

| Objective 2: MUs
We used multiple analyses to evaluate MUs including population structure analyses as above. We evaluated functional connectivity in the Jarbidge-Independence among genetic sites (n = 24) and un-  . Movements were recorded as Euclidean distance of individually marked frogs and could span several years. Frogs in the study area live up to 14 years, but the majority live about 7 years on average (Cayuela et al., 2021). Occupancy data across 1355 sites in the Jarbidge-Independence are derived from surveys conducted between 2001 and 2020 (Nevada Department of Wildlife, unpublished data). Finally, we evaluated hierarchical structure in the neutral data (pruned to markers with <5% missingness) using AMOVA in poppr (Kamvar et al., 2015) across current MUs, watershed units, admixture assignments and assignments based on combined results (e.g., admixture, PCA and functional connectivity).

| Objective 3: Site evaluations and functional connectivity:
We calculated H O , H E and nucleotide diversity from neutral markers using stacks. We estimated effective population size (N e ) for sites with sample sizes ≥11 using the LD method of neestimator version 2.01 (Do et al., 2014). For N e , we pruned site-level data to include neutral, polymorphic SNPs genotyped in at least 80% of individuals, excluding alleles with frequencies less than two critical values, 0.05 and 0.10.
We corrected N e and jackknife confidence intervals based on 13 chromosome pairs in R. pretiosa (Haertel et al., 1974) to reduce bias in N e estimates due to retaining linked loci, which can have a larger effect on N e estimates than either population size or sample size (Waples et al., 2016). We evaluated adaptive differentiation as in Objective 1.
We investigated landscape factors influencing functional connectivity in the Jarbidge-Independence using gravity models (Fotheringham & O'Kelly, 1989;Murphy et al., 2010) implemented in genetit version 0.1-3 (Evans & Murphy, 2020). Gravity models have three components: a measure of spatial proximity (i.e., distance, w), characteristics of the sampled sites that could produce migrants contributing to gene flow (at-site characteristics, v) and characteristics of the landscape matrix restricting or facilitating gene flow between sites (between-site characteristics, c; Table 1). We used (1 -Nei's D) as our measure of gene flow, calculated from the neutral SNPs using adegenet version 2.1.3 (Jombart & Ahmed, 2011). To estimate a singly constrained gravity model accounting for nonindependence of the pairwise genetic matrix, we took the natural log of all dependent and independent variables to linearize the equation and applied a mixed effects model (Murphy et al., 2010). We treated sites as random effects (or grouping factors), accounting for nonindependence in the matrix (Robertson et al., 2018). We sampled the landscape matrix between sites by a straight edge between locations; these edges act as a sample of the landscape and do not represent movement paths. To determine if these edges sufficiently captured environmental variation, we tested multiple edge widths (30 m native cell size, 90, 150, 270 and 510 m) and compared pairwise correlations in independent parameter values between edge widths (Watts et al., 2015). Independent parameter values were highly correlated across edge widths (r > 0.85); therefore, we used the native resolution value (30 m) in all subsequent analyses.
We tested 11 hypotheses of functional connectivity (Table 2). We calculated compound topographic index (Moore et al., 1993), heat load index (McCune & Keon, 2002), relative slope position (De Reu et al., 2013) and surface relief ratio (Pike & Wilson, 1971; data from SRTM digital elevation model, USGS) using spatialeco version 1.3-2 (Evans, 2021), and percentage wetland using landscapemetrics version 1.5.3 (Hesselbarth et al., 2019). Temperature and precipitation metrics were calculated from Rehfeldt (2006) and impervious surfaces were calculated from Dewitz (2020). To assess the spatial context of sampled sites within a larger landscape containing additional wetlands, we At-site Topography, temperature Heat load index hli A measure of incident radiation A proxy for topographic exposure; influences water temperature and primary productivity; high levels reduce productivity in desert regions  Topography Relative slope position rsp Position of a point relative to the ridge and valley, based on its elevation Site productivity is higher closer to ridges in desert regions  and high areas may be inaccessible (Murphy et al., 2010) Moisture Percent wetland p_wetland Percent wetland in a 165-m radius around each site More wetland habitat around a site may increase productivity and facilitate movement (Munger et al., 1998; Moisture; spatial proximity Wetland betweeness wetland_ betweenness Wetlands within a network facilitate stepping stone connectivity The spatial arrangement of wetlands influences overall site productivity (Fortuna et al., 2006; Moisture; spatial proximity Wetland centrality wetland_ centrality Connectedness within a larger network may facilitate long-term stability of a wetland The spatial arrangement of wetlands influences overall site productivity (Fortuna et al., 2006; Source populations; spatial proximity Betweeness, centrality and degree for occupied sites ralu_ betweenness, ralu_centrality, ralu_degree Connectedness to occupied sites may increase the number of potential migrants The spatial arrangement of wetlands influences overall site productivity (Fortuna et al., 2006; Both at-and betweensite Moisture Compound topographic index cti Measure of wetness based on topography (i.e., upslope moisture contribution and ability to hold moisture) Consistently wet at-site conditions have greater productivity; wetter intervening habitat should facilitate connectivity (Munger et al., 1998; Between  (Diestel, 2005), betweenness centrality (Brandes, 2001;Freeman, 1978) and alpha centrality (Bonacich & Lloyd, 2001) for each node for each network. For all continuous variables, we calculated the median value along the edge for each of the landscape variables (e.g., Watts et al., 2015). We tested for correlated variables using a threshold of |r| < 0.70.
Previous work has shown that R. luteiventris exhibits strong isolation-by-distance (IBD; Funk, Blouin, et al., 2005;Murphy et al., 2010;Robertson et al., 2018) making geographical distance the most appropriate null model (Table 2). In addition, spatial autocorrelation is expected in dispersal-limited species such as R. luteiventris and is a condition of the gravity model; therefore, IBD ("distance") is included in all hypotheses. We tested pruned networks at 50 km (the minimum distance to connect the graph across all sites), 100 km maximum distance and saturated (all edges) between observations. Because model results were similar, we present 50-km results as they are more biologically relevant. We used the corrected Akaike information criterion (AICc) for comparison of models fit using maximum likelihood (Akaike, 1973;Burnham & Anderson, 2002). After identifying top model(s), we used restricted maximum likelihood for parameter estimates; if the confidence interval excluded zero, that parameter was considered an important driver of the model (Zuur et al., 2009).
To assess functional connectivity including occupied but unsampled sites (n = 170), we predicted the final model (gravity.predict, population level) to all wetlands with at least one documented observation between 2001 and 2020 (n = 194). To collect covariate data included in the final model, we used the methods above. We then predicted flow at the population level to a 15-km maximum distance graph and a saturated (all edges) graph. We back-transformed the natural log flow predictions using the Naihua correction (Duan, 1983).

| Objective 1: ESUs
We identified 689 candidate adaptive markers on the first two pRDA axes. Most detections were related to precipitation (240 and 222 candidates most strongly correlated with summer and autumn precipitation, respectively), with 227 candidates most strongly TA B L E 2 Eleven hypotheses of Rana luteiventris functional connectivity in the Jarbidge-Independence region and their support based on the Akaike Information Criterion corrected for small sample sizes (AICc), Bayesian Information Criterion (BIC), log likelihood and delta AICc (ΔAICc showing differentiation related to high summer precipitation, and Marys River sites related to low winter temperatures. (Figure 2; Figure S4). WICR in the South Fork Owyhee showed a weaker adaptive signature related to high autumn precipitation and high winter minimum temperatures ( Figure S4). PCA of the 689 candidate adaptive markers did not reveal adaptive divergence associated with the Ruby Mountains until PC7, explaining 2.0% of the genetic variance  Table S3). Population structure analyses identified these regions as distinct, although they were not the primary drivers of genetic divergence among Nevada sites; that is, they did not drive differentiation at the lowest values of K in admixture (Figure 3) or the primary PC axis in PCA (Figure 4).

| Objective 2: MUs
The "best" values of K for admixture and the PCA screeplot were not definitive in the neutral data set ( Figure S3 and S5), due in large part to IBD in the Jarbidge-Independence region (Figures 3 and 4). admixture assignments for K = 9, 11 and 13 illustrate the complex relationships among watersheds ( Figure 5; Figure S6, results for K = 8-13 in Figures S7 and S8 (Table S3), small effective sizes (below) mean that divergence due to drift will require a proportionally larger number of migrants to counteract, in comparison with populations with larger effective sizes (Mills & Allendorf, 1996).
The best-supported model of functional connectivity was the temperature-moisture model ( Table 2,  The genetic site groupings with strongest AMOVA support were delineation using HUC10 watersheds, followed by two assignments based on integration of all results (i.e., admixture, PCA, flow predictions, Table 3; Table S4). These assignments improved among-group variance explained by ~5% over the current watershed-based MUs.
Although HUC10 delineation explained the most among-group variance, this assignment produced 19 groups for 31 sites, creating many singleton sites ( Figure S9). By comparison, assignment into 12 groups based on combined results (Figure 7; Figure S10) explained only 0.1% less among-group variance. The three combined assignments that had 11 groups (Table 3; Table S4) used the 12-group combined assignment (Figure 7) with the following changes: combining Salmon and Marys watersheds ("11 MUs-Salmon+Marys"), no new Winters Creek MU ("11 MUs-SF Owyhee") and no split of the North Fork MU ("11 MUs-North Fork").

| Objective 3: Site evaluations and functional connectivity
Effective population sizes were uniformly small, with an average point estimate of 30.2 (median 22.8, range 3.0-84.7, Table 4; The temperature-moisture connectivity model included the following parameters with 95% confidence intervals that did not overlap zero (Table 6) (Table 2), but the correlation structure of the parameters prevented convergence (producing a singularity error), making the AICc unreliable. In the global model, both connectivity of occupied sites (betweenness_ralu) and wetlands (centrality_wet) had parameter estimates with 95% confidence intervals not overlapping zero.

| DISCUSS ION
Genomic data can play an important role in delineating conservation units, informing our understanding of both neutral and adaptive differentiation in at-risk species. However, these data will rarely be definitive on their own due in part to the difficulty of linking neutral differentiation with demographic independence (Palsbøll et al., 2007;Taylor & Dizon, 1999)

| Evolutionarily significant units within Nevada
The geographical isolation of the southernmost populations of R. luteiventris in Nevada has led to questions about their potential status as ESUs. The Ruby region is 80 km south of the Jarbidge-Independence, with a lack of suitable intervening habitat or hydrological connectivity. The Toiyabe is even more distant, located ~320 km southeast of the Ruby Mountains.
Our genomic results clearly support reproductive isolation of the Ruby and Toiyabe, the first criterion for ESU status, based on neutral pairwise genetic divergence and population structure results.
While population structure analyses do not identify these regions as the primary drivers of genetic differentiation within Nevada, current reproductive isolation is a certainty, given the lack of both intervening suitable habitat and hydrological connectivity. In support of mitochondrial results from a previous study (Funk et al., 2008

| Management units
Management units are generally defined based on demographic independence (i.e., when dispersal rates are <10%; Moritz, 1999;Palsbøll et al., 2007). Translating genetic metrics of population differentiation into estimates of current rates of dispersal (the metric of interest for MUs) is challenging for a number of reasons, in particular because natural systems strongly deviate from many of the assumptions of population genetic models (e.g., Slatkin, 2005;Whitlock & McCauley, 1999). Additionally, methods for inferring contemporary To address these limitations, we combined our regionwide genetic sampling design with occupancy data that included 170 unsampled but occupied sites to project functional connectivity across the Jarbidge-Independence. Projection of the model F I G U R E 5 admixture assignment plots using 39,225 neutral SNPs in 357 individuals distributed across 31 sampling sites. Colours match primary watersheds in Figure 1, except for light green in K = 9 and 13, white in K = 11 and 13, and light grey in K = 11, which indicate genetic differentiation not related to watershed delineation.
to include unsampled sites was essential because Nevada R. luteiventris show limited dispersal capacity, and restriction of the gravity models to the sampled locations would have provided little information on functional connectivity (e.g., only 6% of pairwise Euclidean distances among sampled genetic sites are below the maximum dispersal threshold of 15 km). While this approach does not provide a quantitative evaluation of demographic independence, our evaluation of functional connectivity was sufficient to meet the management goals for the species, which are not highly dependent on an assessment of demographic independence (Taylor & Dizon, 1999). Instead, the goal for MU delineation in Nevada R. luteiventris is to inform recovery efforts that will maintain genetic distinctiveness while maximizing the effectiveness of actions including genetic rescue and reintroductions (Columbia Spotted Frog Technical Team, 2015).
The use of watersheds as a basis for MUs in Nevada was based on previous work that identified ridgelines as barriers to movement for R. luteiventris in Idaho and Montana (Funk, Blouin, et al., 2005;Murphy et al., 2010). However, our connectivity results indicate that movement among Nevada populations may not be highly limited by

| Site evaluations and functional connectivity
The uniformly low effective sizes quantified in our study system are in line with estimates for other western ranids, including northern Fork Owyhee showed a weaker adaptive signature related to high autumn precipitation and higher winter minimum temperatures.
These findings are preliminary without further validation; however, the limited adaptive divergence found across the range indicates that the risk of outbreeding depression across MUs should be minimal during management activities such as genetic rescue Frankham et al., 2011).
Functional connectivity was mediated by temperature and moisture, with higher temperatures reducing at-site productivity and between-site connectivity and wetter conditions increasing at-site productivity and between-site connectivity. This agrees with a recent analysis of 220 sites across the Great Basin, which reported that R. luteiventris occupancy increased with increasing precipitation and decreased with increasing temperature at the watershed scale (Smith & Goldberg, 2020). Support for the importance of connectivity of occupied sites and wetlands in the global model suggests that spatial arrangement in relation to other wetlands, particularly  (Robertson et al., 2018). In contrast, studies of higher latitude populations of R. luteiventris in mountain environments have identified ridgelines and elevation as restricting gene flow (Funk, Blouin, et al., 2005;Murphy et al., 2010) and higher temperatures as increasing at-site productivity (Murphy et al., 2010). These contrasting results reflect how diverse ecological settings alter the landscape template for functional connectivity and productivity, even within a species range (Bull et al., 2011;Robertson et al., 2018).
A core goal of the most recent Conservation Agreement  Figure 1; dual-colour squares and circles represent genomic sites that were reassigned to a new/different MU.
to movement for R. luteiventris in this stronghold of genetic diversity for Nevada populations. This also provides an opportunity to test connectivity models developed using circuit theory, including those projected into the future under different climate scenarios . Predicted flow across the South Fork Owyhee and Independence basins, as well as across the Bruneau and Marys basins represent hypotheses to be investigated in the field (e.g., through habitat evaluation and occupancy searches) as well as in any follow-up genetic work.

| Summary of conservation recommendations for Nevada R. luteiventris
The Jarbidge-Independence region represents the only known reservoir of relatively high genetic variation remaining for Nevada populations. Management actions that maintain a large, functional network of populations and high-quality habitat across this core region could reduce risks to the DPS from stochastic events and environmental fluctuations, including extended drought periods. This will also help maintain genetic diversity and evolutionary potential, which will be important for translocation and genetic rescue efforts to peripheral, genetically depauperate populations.
Due to their smaller N e , lower genetic diversity and isolation, peripheral populations are at much higher risk of reduced fitness due to inbreeding depression and loss of evolutionary potential (Forester et al., in press;Gilpin & Soulé, 1986). These populations are also at high risk from ecological processes such as demographic stochasticity, environmental extremes (particularly long-duration drought events; Pilliod et al., 2021), invasive predatory species and disease . The most serious risks to R. luteiventris populations in peripheral areas differ depending on the basin, but two main concerns are bullfrog invasion and habitat loss/degradation due to anthropogenic alterations to regional hydrology and climate TA B L E 4 Effective population size (N e ) estimates using a maximum allele frequency of 0.05 (i.e., Pcrit = .05). N e and jackknife confidence intervals (CIs) are corrected for the number of chromosome pairs in Rana pretiosa. The "inf" estimate indicates that the confidence interval includes infinity. change. In locations where bullfrogs are thought to be causing population declines and extirpations, managers could consider translocations to bullfrog-free habitats as an effective, primary option. Where possible, translocations that include genetic rescue could also enhance genetic diversity and evolutionary potential (Bell et al., 2019;Ralls et al., 2018Ralls et al., , 2020, work that is currently in progress in the western part of the Jarbidge-Independence based on these results.
In regions where maintaining genetic distinctiveness is important, selecting source populations based on genetic similarity (e.g., from PCA, admixture and predicted flow results) and/or geographical proximity would be helpful, even if those populations have reduced genetic variability. However, where maximizing genetic rescue and fitness effects is most important, selecting sites with higher diversity (i.e., from the core of the Jarbidge-Independence) would be best, especially because risks of outbreeding depression should be low (Bell et al., 2019;Frankham, 2015).
Where habitat loss and degradation are the driving factors in population declines, improvements in habitat quality that restore natural hydrology can clearly boost population sizes, as shown by restoration efforts in the Toiyabe region and elsewhere Pilliod & Scherer, 2015). Climate change in the Great Basin is already increasing winter temperatures and aridity, and expanding the frequency and duration of droughts (Snyder et al., 2019).
Restoration of resilient hydrological regimes will be incredibly diffi- population sizes, it will not restore genetic diversity in the absence of genetic rescue efforts, as discussed above.

| General recommendations
Our study highlights both the strengths and the limitations of genomic data in informing diverse management questions in complex natural systems, and how integration of genomic with complementary data can improve downstream inferences. First, by carefully designing our genetic sampling across the Nevada range, we were able to identify unexpected patterns of population connectivity across watershed boundaries that highlight the uniqueness of desert populations of this widely distributed species.
This genetic sampling design required the unavoidable trade-off between sampling more sites and genotyping more individuals per site. In our case, we sacrificed sample size at any given site to improve geographical and environmental coverage across the range, a design that allowed us to better meet the diverse management goals of the study. Despite careful planning, we still have areas of uncertainty that will require follow-up if further resolution is needed (e.g., unassigned units in Figure 7). Transparency about uncertainty and the limitations of sampling designs and analytical inferences is an essential component of both effective communication with management partners and application of study results to conservation actions.
Second, we took advantage of complementary data that included capture-recapture and range-wide survey data to improve our inference of functional connectivity and MU delineation in the absence of information on demographic connectivity. A dispersal threshold derived from capture-recapture data allowed us to apply a biologically relevant limit on the connectivity model, improving its relevance for delineating MUs and highlighting movement corridors that are more likely to be ecologically relevant. Similarly, site occupancy data allowed us to predict the gravity model beyond the small subset of genomic sampling sites, expanding our capacity to infer connectivity across the landscape. This approach provides a framework for MU delineation in other species where funding limitations either prevent genotyping of all occupied sites or do not allow for sufficient sample sizes to evaluate contemporary dispersal using genomic analyses. More generally, communication with managers on the functional definition of an MU that best meets the species management goals is critical to designing an effective and beneficial conservation genomics study.
Finally, genomic-scale data allowed us to identify adaptive differentiation, or in this case, the relative lack of it, which is useful for informing ESU delineation and genetic rescue efforts. A lack of genomic resources in many at-risk species, such as high-quality, annotated reference genomes, combined with the practical difficulties of experimental approaches for demonstrating local adaptation in fragile, declining populations (especially in cases where generation times are long) will continue to hamper the utility of genomic assessments of evolutionary potential in at-risk species .
Despite these limitations, anonymous genomic assessments of candidate adaptive differentiation, such as the approach used here, remain useful in identifying spatial patterns of differentiation to inform CU delineation as well as reducing the risk of potential outbreeding depression during genetic rescue. Overall, while a single genetic data set is unlikely to be able to fully address all potential conservation and management questions for a given species, a carefully planned study, designed in collaboration with practitioners to address priority management questions, can provide diverse insights for conservation and recovery actions.