Broad‐ and fine‐scale structure across the distribution of the relict dace (Relictus solitarius) in the Great Basin desert, USA

The Relict Dace is a small cyprinid and the only species in its genus, Relictus. It is naturally distributed in four drainage basins of central Nevada ‐ Butte Valley, Ruby Valley, Steptoe Valley, and Goshute Valley. The species has experienced periods of isolation and connection within these four basins since the late Pliocene, with the last 100 years characterized by anthropogenic disturbance. To better inform conservation actions, we investigated range‐wide and intra‐basin genetic structure and diversity and conducted Fst outlier tests using RAD‐sequencing. We found high levels of genetic structure and four main lineages on a broad scale corresponding primarily to geography, and on a very fine scale (<1 km) within drainage basins. Signatures of selection in the form of outlier loci were documented in multiple locations across the range. Our results provide baseline data for conservation efforts and highlight the degree of fine‐scale genetic structuring that may be present in fishes of the Great Basin that are now restricted to small, isolated habitats—particularly wetland complexes. We caution against assuming that geographic distance is a proxy for genetic similarity in Great Basin fishes, encourage thorough genetic sampling, and note that broad‐scale genetic surveys can overlook fine‐scale genetic structure.


| INTRODUCTION
Many aquatic species naturally occur in isolated and infrequently connected populations unevenly distributed over a landscape. The persistence of such populations can reach anywhere from shallow to more ancient time scales, with many thousands of years of isolation among populations of the same species. Species or species complexes that are characterized by multiple isolated populations often lie at the intersection of both evolutionary interest and conservation concern, particularly in arid landscapes such as the Great Basin Desert (Meffe & Vrijenhoek, 1988;Sa glam et al., 2016;Smith, 1978). Freshwater fishes in the Great Basin present numerous cases of vulnerable isolated populations that exhibit physiological adaptations to extreme conditions such as the Devils Hole Pupfish (Cyprinodon diabolis), where the entire species is restricted to a single water-filled cavern with constant 33 C temperatures (Sa glam et al., 2016). More recently, human influences including groundwater pumping and introduced species have broadly threatened native aquatic fauna in the Great Basin (e.g., Minckley & Deacon, 1991;Unmack & Minckley, 2008). The history of many Great Basin fish species has been one of survival and adaptation together with extinction and extirpation as they persist through increasing perturbations to their unique and limited habitats (Minckley & Deacon, 1991).
The Relict Dace (Relictus solitarius) is one such unique extant Great Basin fish and is the only species in the genus Relictus. It is endemic to and the only fish naturally found in four adjacent closed basins in the northeastern portion of the state of Nevada (Goshute, Steptoe, Ruby, and Butte Valleys) and was introduced into a fifth (Spring Valley). Although Pleistocene pluvial conditions connected much of neighboring valleys and fish assemblages, geologic, and electrophoretic evidence shows that Relict Dace evolved pre-Pliocene and remained isolated in its four native drainages throughout the Pleistocene (Hubbs & Miller, 1948). Unlike many other Great Basin Fishes which inhabit lakes and river systems, most Relict Dace populations occur in spring systems that originate from a very limited number of springheads with highly restricted outflow channel habitats. Hubbs and Miller (1948) described Relict Dace as "rather degenerate" relative to the closely related Rhinichthys, attributing these traits to the lack of swift water habitat and aquatic predators. These characteristics also make Relict Dace vulnerable to habitat alteration, predation, and competition from introduced fish; numerous extirpations and threats have been recorded since at least the 1980s (Hubbs et al., 1974;Vigg, 1982).
Throughout their range, Relict Dace show very high levels of structure and diversity and more recent admixture in the eastern valleys when compared to the western valleys (Houston et al., 2015), likely due to more available habitat in Goshute and Steptoe Valleys. Previous mtDNA data suggest an initial split between the eastern (Goshute and Steptoe) and western (Butte and Ruby) Valleys $1 million years ago, followed by significant diversification during the Pleistocene and increased isolation through the Holocene from drying, and more recently through habitat alteration and introduced species (Houston et al., 2015). These findings and observations of other spring obligate species, such as springsnails (Hershler et al., 2014) also suggested unsampled diversity and connectivity patterns in Relict Dace, particularly in the larger wetland systems such as Johnson Springs Wetland Complex (JSWC) in the eastern valleys. Though many Relict Dace populations are currently disconnected-even within single spring systems-the large wetland complex of the JSWC is a notable exception. JSWC is a uniquely large, ecologically diverse, and largely intact wetland complex with no introduced fish species located at the northern end of Goshute Valley (Figure 1). While most Relict Dace live in very small single-spring locations, in contrast JSWC is approximately 2.4 km in length (north-south) with a width of 1.1 km and a perimeter of approximately 6.3 km. A total of 88 discrete springs have been identified at JSWC, and these spring types have further been classified as hillslope, helocrene, limnocrene, and hypocrene habitats (SWCA Environmental Consultants, 2020). Though in 1982 it was believed that Relict Dace at Big Springs Ranch (JSWC) in Goshute valley were extirpated (Vigg, 1982), Relict Dace currently occupy 29 separate springs in addition to a number of man-made ditches and channels that were historically (within the last $100 years) used for irrigation purposes in addition to an ephemeral reservoir on the property. The size, diversity of habitat types, amount of available habitat, and resulting densities of Relict Dace at JSWC is greater than at any other existing location, in many cases by an order of magnitude (Crookshanks, 2006;SWCA Environmental Consultants, 2020). Combined with the notable lack of nonnative species which threaten populations in the other valleys, JSWC is considered the most secure habitat for Relict Dace across its entire range (Crookshanks, 2006).
At the time of range-wide sampling for Houston et al. (2015), managers were primarily concerned with hybridization between Speckled Dace and Relict Dace in Butte and Ruby valleys, and it was assumed that most of the JSWC were occupied by a single population which could be represented by Big Springs. However the high habitat diversity of JSWC, results on the range-wide distribution of genetic variation from Houston et al. (2015), and recent comprehensive surveys uncovering significant habitat diversity prompted managers to question (1) Is JSWC is a single large population or multiple populations, (2) Is there evidence of local adaptation in the form of outlier loci? and (3) to What extent are populations within JSWC connected, if at all?
In this paper, we characterize evolutionary relationships across the range of Relict Dace and examine population genetic relationships within the major geographic divisions of the species. We increase both genetic and geographic sampling over previous studies, refining the broad strokes of Relict Dace relationships and providing hypotheses of genetic structuring at a fine scale. We leverage next-generation sequencing to examine how both neutral and adaptive genetic variation is distributed in Relict Dace throughout their range. In particular, we focus on JSWC to examine population genetic structure at a fine-grained spatial scale and explore evolutionary processes that may contribute to genetic subdivision within JSWC.  Figure 1). Locations used in Houston et al. (2015) are noted, as well as alternative names for locations. Fin clips were dried and stored in envelopes until DNA extraction using the Qiagen DNeasy kit (Qiagen) according to the manufacturer's protocol.
Five RAD libraries were prepared following the protocol described in Ali et al. (2016) using 150-200 ng of DNA from each individual and the restriction enzyme Sbf1. Libraries were pooled into a single lane for paired-end 150 base pair (bp) sequencing on an Illumina HiSeq4000.

| De novo reference assembly
RAD sequence de-multiplexing and de novo reference assembly were conducted according to the pipelines outlined by Miller et al. (2012) andSa glam et al. (2016). We selected seven Relict Dace from Franklin Lake (Ruby Valley) with the greatest number of raw reads (1,811,832-2,353,000 reads) to generate a de novo reference contig assembly, using the resulting alignment to identify 31,390 putative RAD loci, which were then extended into longer contigs using paired-end data and the genome assembler PRICE (Ruby et al., 2013).  Next, we filtered contigs for a minimum length of 300 bp, retaining a set of 28,324 RAD loci ranging in size from 300 to 738 bp. The final set of RAD loci was combined into one FASTA file, indexed using the Burrows-Wheeler Alignment (BWA) tool , and used as a de novo partial genome reference for all subsequent analyses.

| Alignment and paralog filtration
De-multiplexed sequences for all individuals were aligned to our de novo set of reference RAD contigs using the BWA-minimal exact match (MEM) algorithm implemented in the BWA software . Alignment maps were outputted as Sequence Alignment/ Map (SAM) files and converted to Binary Alignment/ Map (BAM) files using SAMTOOLS (Li, 2011;. We used SAMTOOLS to sort, filter for proper pairs, remove PCR duplicates, and index BAM files. Individuals with fewer than 100,000 aligned reads after all filters were removed from further analysis. We used the program NGSPARALOGS (https://github.com/tplinderoth/ngsParalog) to identify and remove putative paralogous loci.

| Population genetic structure and phylogenetic relationships
We examined population genetic structure with three different sample sets of Relict Dace, (1) range-wide (all five valleys), (2) western valleys (Butte and Ruby), and (3) eastern valleys (Steptoe, Spring, and Goshute). All population genetic analyses were conducted using the probabilistic framework implemented in ANGSD (Korneliussen et al., 2014). Analyses were conducted using the GATK genotype likelihood model (GL 2) with a minimum mapping quality of 10 (minMapQ) and a minimum base quality of 20 (minQ). We identified polymorphic sites (SNP_pval 1e-6), inferred major and minor alleles (doMajorMinor), estimated allele frequencies (doMaf 2), and retained sites that were present in at least 50% of individuals and had a minor allele frequency of at least 0.05 (minMaf). As the presence of closely-related individuals may lead to false signatures of population genetic structure (reference), a kinship matrix was calculated with PCANGSD (Meisner & Albrechtsen, 2018) (kinship). We then performed principal component (PC) and admixture analyses with PCANGSD using the kinship matrices masking individuals with a coefficient of relatedness (r) ≥0.0625 (third degree), the default setting of PCANGSD. Principal component and admixture results were visualized using the tidyverse collection of R packages. Locations sampled by Houston et al. (2015).
We constructed a range-wide phylogeny by utilizing the range-wide sample list with closely related individuals removed as indicated by PCANGSD. Genotype calls were generated with ANGSD through the creation of a plink-formatted file (doPlink 2) while retaining sites in 95% of individuals (minInd 344), restricting minor allele frequency to be at least 0.05 (minMaf 0.05), and a SAMtools genotype likelihood model (GL 1) with a posterior cutoff of 0.95 (postCutoff 0.95). We inferred major and minor alleles from genotype likelihoods (doMa-jorMinor 1) and calculated major and minor alleles with a fixed major and minor (doMaf 1) using frequency as a prior (doPost 1). A minimum mapping value of 20 (minMapQ 20) and minimum base quality of 20 (minQ) were also specified. The resulting genotype calls were re-coded to a variant call format (VCF) file with PLINK (Purcell et al., 2007). We then removed linked SNPs with the BCFtools prune with the following specifications: bcftools +prune -l 0.9 -w 10,000 (Danecek et al., 2021). The pruned VCF file was converted to a PHYLIP file with vcf2phylip.py (https://github.com/ edgardomortiz/vcf2phylip) and ascertainment bias correction performed with ascbias.py (https://github.com/ btmartin721/raxml_ascbias/blob/master/ascbias.py). Individuals with missing data comprising more than 10% of the sequences were removed with removeMissing.pl (in this manuscript's code repository). A species tree of sampling locations was constructed with SVDQuartets within PAUP* by evaluating 1,000,000 random quartets and support assessed with 100 bootstraps (Chifman & Kubatko, 2014Swofford, 2003). The tree was visualized with the ggtree package in R, and node supports <70 were collapsed within R using the as.polytomy function of ggtree (Yu et al., 2017).

| F and neutrality statistics
We investigated genetic differentiation between populations and within groupings ascertained with previous bioinformatic steps by calculating pairwise F ST values using all individuals. To do this, we first estimated an unfolded site frequency spectrum (SFS; doSaf) for each population and then used REALSFS (Korneliussen et al., 2014) to calculate the two-dimensional SFS and global F ST value (weighted) for each pair of populations.
To measure relative levels of genetic diversity within individual Relict Dace collection sites, we calculated the average number of pairwise nucleotide differences (nucleotide diversity; ϴπ) and the number of segregating sites (Watterson's theta, ϴ W ), including related individuals. First, a folded SFS was estimated for each population and used as a prior (pest) to determine per-site thetas (doThetas). These values were then averaged across all loci to produce ϴπ and ϴ W estimates. We then computed Tajima's D, the difference between ϴπ and ϴ W scaled such that they would be equal in a neutrally evolving population of constant size. A negative Tajima's D could signify a population expansion or purifying selection, and a positive Tajima's D could suggest population contraction or balancing selection.

| Outlier loci
To identify outlier loci in populations while accounting for drift due to population structure and deviations from Hardy-Weinberg due to inbreeding, we conducted genomewide selection scans on all individuals based on population structure using PCAngsd (selection 1). Based on the FastPCA model (Galinsky et al., 2016), PCAngsd performs a selection scan along all significant PCs and identifies variants (SNPs) whose differentiation along the top PCs is significantly greater than the null distribution expected from genetic drift. It does this by minimizing confounding effects due to inbreeding by filtering out SNPs that show significant deviations from HWE (inbreedSites) based on the kinship matrix between all individuals (inbreed 3). The resulting selection statistic is χ 2 -distributed with 1 degree of freedom. This statistic is then used to identify outlier variants along each significant PC using a Bonferroni correction (based on the number of significant PCs and sites). Outlier loci are those regions that show distinct deviations from Hardy-Weinberg due to selection and not drift, though we cannot assess whether these outliers are from the past or current selection, or rule out false negatives.

| RESULTS
Sequencing yielded a total of 499,659,476 raw reads, amounting to an average of 1,209,368 raw reads for 396 Relict Dace individuals. From the de novo assembly, we produced a de novo reference consisting of 28,324 RAD loci with an average length of 421 bp. The average alignment success to this reference was 89.1%. See Table S1 for individual alignment statistics. Four individuals with <100,000 filtered aligned reads were removed from further analysis.
F I G U R E 2 Range-wide genetic relationships of relict dace. For the range-wide data set, relationships among individuals are shown through principal component (PC) analysis, with the first two PCs shown in 2A and the first and third in 2B with the percentage of the variance explained for each PC indicated on the axis label. The admixture plot of K =2 genetic clusters is shown in 2C and K =3 in 2D with the individuals placed along the x-axis with sampling locations indicated and admixture coefficient (Q) on the y-axes; K =2 and K =3 highlight the early division between the eastern and western valleys and Twin Springs. The optimal K = 16 value for all populations together is shown in Figure S1. A species tree hypothesis for relict dace with sampling locations indicated at tips with point size proportional to the number of individuals included in the phylogenetic analysis. Values at nodes indicate bootstrap support values, and those <70 were collapsed. In all panels, samples are color-coded by geography; red = western valleys, blue = eastern valleys excluding Twin Springs, green = Twin Springs. Further detail for sampling locations is available in Table 1 7, RBRL_246 7, and JSWC BS 4, CS04 2, CS16 2, CS22A 1, and SS19 6. Considering all (n = 392) Relict Dace samples analyzed, an initial 84,591 SNPs were assigned genotype likelihoods for the range-wide data set. After the removal of the 30 closely-related fish and repeated MAF filtering, 82,936 SNPs were available for PC and admixture analyses. For phylogenetic analysis, an initial 30,673 called SNP genotypes were pruned to 26,700 SNPs. Removal of individuals with substantial missing data (>10%) produced a data set of 338 individuals and 21,395 SNPs after correcting for ascertainment bias. The input data file and output tree file are available at https://github.com/MacCampbell/relictdace/blob/master/svdQuartets-files/.
Overall, PC, admixture, and phylogenetic analyses of the range-wide dataset divided Relict Dace into two distinct groups: western valleys (Butte and Ruby) and eastern valleys (Goshute, Steptoe, and Spring) (Figure 2). The first PC axis (representing 17.38% of the variation) separates western valleys from the eastern valleys (Figure 2a), and PC 2 (4.62% of the variation) separates individuals from Twin Springs in Goshute Valley (GSTS) separated from all other Relict Dace in the eastern valleys. Admixture analysis produced similar results (Figure 2b). At K =2, the two observed clusters also split the western and eastern valleys, and at K =3 GSTS separated from the eastern valleys to form the third cluster (Figure 2c,d).
The optimal K as determined by PCAngsd is 16 and is presented in Figure S1. Phylogenetic analysis indicates Butte Valley locations are the earliest-branching followed by the Ruby Valley locations (Figure 2e). There is maximal (bootstrap = 100%) support for the monophyly of the eastern group (Goshute, Steptoe, and Spring). Within the eastern valleys, GSTS is placed as the earliest-branching lineage with monophyly of JSWC in Goshute Valley receiving maximal support. Steptoe Valley and Spring Valley receive low support for monophyly (bootstrap = 75%) and Spring Valley is nested within Steptoe Valley.

| Western and Eastern Valley substructure
To examine substructure within the western and eastern valley lineages individually, we again used PC and admixture analyses, excluding GSTS from the eastern valleys due to its distinctiveness. The western valleys again were examined for related individuals and a single individual was removed from BTQP, reducing the 100 individuals to 99. The 50,633 genotype likelihoods initially generated were reduced to 50,345 after the removal of the single BTQP individual and the application of MAF filtering. In the western valleys we found four clusters of Relict Dace: (1) Odger's Creek (BTOC) and Quilici Pond (BTQP) (2) Franklin River (RBFR) and Franklin Lake (RBFL) (3) Ruby Lake Pond 226 (RBRL_P226) and (4) Ruby Lake Pond 246 (RBRL_P246) (Figure 3a-c). See Figure S2 for plots of all K values examined in the western valleys.
The number of genotype likelihoods generated from the eastern valleys was, including all individuals, 87,712 from 257 individuals. Filtering out 14 closely-related individuals and again with MAF filtering reduced this number to 85,931. Individuals were removed from: BS 4, CS04 2, CS16 2, CS22A 1, and SS19 5. In the eastern valleys, we found a substantial population substructure (Figure 4). Steptoe Valley, genetic differentiation between sites was high (F ST =0.08-0.23), and STMR was the most differentiated.
With the exception of GSTS, diversity was higher in the populations in the eastern valleys (mean ϴπ = 0.0030-0.0066; Figure 3e) than the western valley populations (mean ϴπ = 0.0023-0.0034; Figure 4e). GSTS had slightly lower diversity levels than other eastern valley locations (ϴπ = 0.0030) (Figure 4e). The highest ϴπ values were found in the JSWC admixed locations NS03, NS05, and SS19, and the large STRT, and STWA locations (Figure 4e). Like ϴπ, ϴ W values in the western valleys (mean ϴ W =0.0021-0.0031) were lower than in the eastern valleys (mean ϴ W =0.0025-.0069). In the eastern valleys, GSTS had the lowest value (0.0025) and the JSWC admixed locations plus STRT and STWA had the highest values ( Figure S4). Tajima's D values ranged from À0.19 to 0.68 in both the eastern and western valleys, suggesting that none of the Relict Dace populations have recently experienced large deviations from neutrality, or detectable large demographic changes (i.e., population expansions or bottlenecks). All diversity statistics are presented in Figure S4.

| Outlier loci
Based on the major division uncovered between the eastern and western valleys, we ran selection scans along all significant PC axes obtained from two separate PCA results. For the western valleys, we interrogated three PC axes and 48,444 SNPs resulting in a genome-wide significance threshold of 3.44 Â 10 À7 (α =0.05; 48,444 SNPs Â 3 PCs). For the eastern valleys, we interrogated 13 PCA axes and 84,240 SNPs resulting in a genome-wide significance threshold of 4.56 Â 10 À8 (α = 0.05; 84,240 SNPs Â 13 PCs).
In the western valleys, selection scans conducted along the three significant PC axes captured eight outlier sites on PC 2 and 4 sites on PC 3 (along 12 RAD loci) that reached genome-wide significance (LTR ≥25; p-value <3.44 Â 10 À7 ) ( Table S2A) (Table S2A).

| DISCUSSION
The Great Basin is known for its isolated populations of aquatic fauna that have experienced reticulation events throughout the Pleistocene and drying Holocene (Smith et al., 2002). Relict Dace occupy a unique place in this system, as geologic evidence suggests that unlike many other Great Basin fish taxa they were not connected to adjacent basins and evolved into a monotypic genus, Relictus. Using genome-wide data and intensive surveying, we examined the genetic structure and diversity, as well as potential signals of local adaptation in the form of outlier loci from Relict Dace taken from locations throughout their native and introduced range. Our goal was to provide relevant information for managers intent on preserving biodiversity in Relict Dace. Our analyses covered broad and fine spatial scales, with a particular focus on Johnson Springs Wetland Complex (JSWC). At broad (inter-basin) spatial scales, our expectations were met and we found a divergence between locations to agree with neutral expectations of genetic divergence through physical isolation. However, on a finer scale, we found both high levels of structure and loci putatively under selection in individual populations.

| Range-wide structure
The genetic structure and differentiation between Relict Dace sampling locations are largely hierarchical in nature, with greater divisions at larger (inter-valley) spatial scales compared with small (intra-valley) spatial scales. Our observations are consistent with the evolutionary history of other Great Basin fishes where large, genetically diverse populations were fragmented by the recession of Pleistocene pluvial lakes (Billman et al., 2010;Smith et al., 2002).
The results presented here demonstrate this hierarchy and complex history with a clear resolution of inter-valley phylogeographic relationships and less well-resolved intravalley relationships (Figure 2). The principal component analysis identifies three major divisions of Relict Dace, broadly western valleys, eastern valleys, and Twin Springs (GSTS). These three genetic clusters are recapitulated by admixture and phylogenetic analysis. The increase in the number of variable and independent sites compared to Houston et al. (2015) results in a concomitant improvement in terms of phylogenetic resolution of eastern valley sampling locations (Figure 2e). Excluding GSTS, the eastern valley sampling locations resolve into a strongly supported monophyletic JSWC lineage (bootstrap = 100%) that is sister to a moderately supported Steptoe Valley lineage (bootstrap = 75%).

| Western valleys
Butte and Ruby valleys are thought to have been hydrologically connected more recently (possibly mid-late Pleistocene) by Hubbs et al. (1974), and this is supported by their genetic similarity. In Ruby Valley, Relict Dace occurs in the Ruby Lake and Franklin River systems. The two locations sampled in the Ruby Lake system were ponds RBRL_P226 and RBRL_P246, which are in the Ruby Lake National Wildlife Refuge where there is active management to remove nonnative species. Yet the ponds are genetically divergent, and we detected outlier loci in each. In the hydrologically connected Franklin Lake and Franklin River, we found the locations to be genetically similar and did not detect outlier loci, though we cannot rule out non-neutral processes. Our Butte Valley samples consist of locations along Odgers Creek and Quilici Pond, which did not display genetic structuring (Figure 3). Quilici Pond is located 107 m northwest of Odgers Creek and is thought to connect during high water events (H. Korell, pers. comm.). The overall genetic similarity of Odgers Creek and Quilici Pond may originate from current or recent gene flow or more recent colonization of one water body to another.

| Twin Springs
The eastern valleys-Goshute, Steptoe, and Spring-are divergent from each other, though less starkly than Ruby and Butte valleys, and our results indicate more recent gene flow. GSTS is located in the southern end of Goshute Valley between JSWC and Steptoe Valley, but our analyses do not indicate it is most closely related to either JSWC or Steptoe Valley Relict Dace (Figure 2e). Previously with mtDNA data, Houston et al. (2015) found most fish from GSTS to have the same haplotype and the sampling location was monophyletic and nested within generally polyphyletic JSWC and Steptoe Valley sampling locations. The complex evolutionary history of GSTS including a loss of diversity due to drift and the rapid separation of Relict Dace into refugia results in few shared sites to inform branching patterns, and the few phylogenetically informative sites that do support a "true" topology may at times be less abundant than the sites that support a "false" topology. Rapid separations of populations or species are difficult to resolve, e.g., Harrington et al., 2021). We find that ambiguity exists in the data regarding the relationships among major eastern valley phylogeographic divisions as indicated by nonmaximal bootstrap support values (Figure 2e). A sister lineage of GSTS within the eastern valleys may be resolved with an additional geographic sampling of nearby locations and applying the multispecies coalescent.

| Johnson Springs wetland complex
In JSWC, our analysis uncovered striking population genetic structure as well as evidence of predominantly one-way gene flow and potential for local adaptation. Our results support seven distinct genetic clusters from 10 sampling locations within JSWC. The seven distinct genetic clusters are restricted to single springheads (BS, CS02, CS03, CS04, CS11, CS16, and CS22A) which roughly "flow" intermittently toward three admixed locations (NS03, NS05, and SS19). Three unique locations (CS03, CS04, and CS11) have loci with locally elevated minor allele frequencies that may suggest the selection. The admixed locations NS03 and NS05 are located on the opposite side of the complex from SS19, but all three appear to be "genetically downstream" of the unique springheads. There are several explanations for one-way gene flow in JSWC. One is that some locations (specifically some of the pothole springheads) would be difficult to return to due to very low or indiscernible outflow (Crookshanks, 2006). Another possible cause is the ephemerality of the locations not directly fed by permanent springheads. It is not uncommon for surveys of relict dace to find no detectable fish in survey sites where they were recently present only to find that fish are abundant in a subsequent year. Importantly, the fish captured in the proximate NS03, NS05 locations and the farther away SS19 have been admixed for at least several generations and show similar admixture proportions; i.e., there is no dominant ancestral proportion. It is possible local extirpations happen infrequently, or fish are able to find refuge in other locations and then return. Stochastic extirpations caused by a variety of factors happen locally in wetland complexes such as JSWC, even in springheads. In CS03, an elk died in 2017 and wiped out the population after the samples analyzed here were collected.
The relative levels of diversity and differentiation that we found in JSWC are expected given admixture results: springhead locations have lower diversity and higher differentiation and the non-springhead locations have higher diversity and are less differentiated. One unexpected finding is low diversity in BS, which is the largest springhead habitat with the highest levels of outflow in JSWC, though BS is sourced from a different aquifer than the other springs at JSWC.
The hydrology of JSWC is spatially and ecologically complex, with variability on a seasonal and annual basis. Hitherto this study, the channels that cross the JSWC were thought to provide ample opportunity for gene flow and admixture, leading biologists to hypothesize low levels of structure. Though our study uncovered unexpected structure and potential for differing selective pressures, much of the complex remains genetically unsampled, specifically channels that connect some of the locations sampled herein. Further surveys in JSWC would shed additional light on connectivity and ecological and phenotypic data on selective processes and provide information relevant to management for one of the most important locations for Relict Dace.

| Spring Valley
Though records are incomplete, the two Spring Valley locations of SPSH and SPKR is thought to be founded from a location in Steptoe Valley. Our analysis strongly supports this hypothesis, as the Spring Valley locations are nested within the Steptoe Valley Clade (Figure 3). Yet we did not identify an obvious source location among our samples. The relatively low diversity estimates for Spring Valley locations are concordant with a founder event (Figure 4e and Figure S5). The uncertainty on the source populations for SPSH and SPKR are likely because either the source was not sampled, or demographic effects have obscured the source.

| Steptoe Valley
The four Relict Dace locations we sampled in Steptoe Valley cover a much broader area compared with JSWC, and for the most part we observed typical Great Basin broadscale intra-valley patterns where populations are most closely related within the valley (Smith et al., 2002;Figure 2e). Our PCA grouped three of four locations (STWA, STRT, and the two ponds at STPR) closely together ( Figure 4a) and these four locations exhibit high diversity levels (Figure 4e), suggesting relatively more recent connections and large population sizes, though none of these groups appears to have very recent gene flow.
Steptoe Valley Wildlife Management Area (STWA), and STRT are both large wetland complexes composed of ponds and ditches connecting them. STWA is one of the original sites sampled by Miller in 1938 (Hubbs &Miller, 1972). STRT is on private land used for agriculture and has abundant exotic species. Though surveys for this study found that many of the STWA ponds had been invaded by Utah Chub (Gila atraria), the two ponds sampled for this study still hold substantial numbers of Relict Dace and levels of genetic diversity.
The most differentiated population in Steptoe Valley is STMR (Morrison Ranch) despite its geographic location between STPR (Phillips Ranch) and STWA/STRT. However, the phylogenetic tree places STMR sister to STPR within the clade composed of Steptoe Valley sampling locations. Both STPR and STMR (Morrison Ranch) are located on private property in relatively large wetland complexes. The two locations sampled at STPR are located $6 km from each other. A close relationship between STPR and STMR is expected because prior to human settlement ($100 ya) waterways connected to STPR and STMR discharged into the Duck Creek drainage which flows from south to north; outflow from STMR would have flowed toward STPR. It is unlikely that there are Relict Dace in Duck Creek today due to the presence of predatory Northern Pike (Esox lucius) and Largemouth Bass (Micropterus salmoides).

| Outlier loci
For Relict Dace and other species inhabiting ecologically diverse habitats, it is not unexpected to find outlier loci. Just within JSWC, periodic surveys have recorded Relict Dace in areas with temperatures ranging from 2.2 to 22.4 C (SWCA Environmental Consultants, 2020). Our identification of loci putatively under selection is congruent with the observation of high phenotypic variability within the Relict Dace across its range (Hubbs et al., 1974). Future research with an annotated genome and phenotypic and environmental data may shed light on the selective agents in play, but these findings should not be used by managers until a greater understanding can be presented.

| CONCLUSION AND CONSERVATION IMPLICATIONS
Anthropogenic landscape alteration of Relict Dace habitats within the last $100 years complicates our understanding of population history and management, as humans have directly and indirectly facilitated and hindered movement across the landscape. Our work highlights the utility of extensive regular genetic and physical surveys of these arid wetland complexes, though it requires careful planning and can be physically and logistically difficult. Many of the locations are situated on private land with active agricultural use where access is denied. Even on public land surveying can be difficult; Relict Dace can live in shallow channelized habitats where emergent vegetation biomass limits open water. Furthermore, wetland complexes are highly temporally variable, to the extent where some locations do not reliably hold detectable Relict Dace on an annual basis. In addition, many of the sampled locations including Steptoe and Ruby Valley National Wildlife Refuges hold introduced species which can result in stochastic population declines or extirpations. Though JSWC is the most secure source of Relict Dace based on its size and lack of introduced fish species, its hydrology is affected by groundwater pumping associated with an adjacent active gold mine. Despite these obstacles, we recommend that to the extent possible biologists continue to survey Relict Dace populations. While many locations with Relict Dace may never be secured and some have been extirpated, the larger national wildlife refuges and JSWC hold a significant amount of diversity and should be protected. performed all laboratory work and conducted initial data analysis and writing. Matthew A. Campbell and Ismail K. Sa glam contributed to data analysis and writing. Chris Crookshanks provided critical review and coordinated site visits and sample collection.