Characterizing patterns of genomic variation in the threatened Utah prairie dog: Implications for conservation and management

Abstract Utah prairie dogs (Cynomys parvidens) are federally threatened due to eradication campaigns, habitat destruction, and outbreaks of plague. Today, Utah prairie dogs exist in small, isolated populations, making them less demographically stable and more susceptible to erosion of genetic variation by genetic drift. We characterized patterns of genetic structure at neutral and putatively adaptive loci in order to evaluate the relative effects of genetic drift and local adaptation on population divergence. We sampled individuals across the Utah prairie dog species range and generated 2955 single nucleotide polymorphisms using double digest restriction site‐associated DNA sequencing. Genetic diversity was lower in low‐elevation sites compared to high‐elevation sites. Population divergence was high among sites and followed an isolation‐by‐distance model. Our results indicate that genetic drift plays a substantial role in the population divergence of the Utah prairie dog, and colonies would likely benefit from translocation of individuals between recovery units, which are characterized by distinct elevations, despite the detection of environmental associations with outlier loci. By understanding the processes that shape genetic structure, better informed decisions can be made with respect to the management of threatened species to ensure that adaptation is not stymied.

| 1037 GIGLIO et aL. translocations, used to bolster declining populations, can be improved by incorporating genetic data to tailor translocation actions to outcomes that boost genetic variation. Selecting source populations that are genetically appropriate for the target population (Johnson et al., 2010), prioritizing target populations with low genetic variation (Whiteley et al., 2015), or gauging incorporation of source genotypes (Bateson et al., 2014;Latch & Rhodes, 2005;Mulder et al., 2017) can improve efforts to retain genetic variation. Conservation and management actions could further improve the evolutionary potential and adaptive capacity of populations by incorporating data from studies of adaptive variation, especially if those actions include translocations, genetic rescue, or assisted gene flow (Flanagan et al., 2018;Funk et al., 2018).
Advances in genomics for nonmodel species means that researchers can generate broad coverage and high-resolution genomic data for an increasing number of species. Genomic data can be used to survey both adaptive and neutral genetic variation and can be incorporated into conservation policy to improve long-term species viability against changing environments and exposure to new diseases (Flanagan et al., 2018;Funk et al., 2018).
In this study, we use a population genomics approach to understand the maintenance of genetic variation in the threatened Utah prairie dog (Cynomys parvidens). The Utah prairie dog is one of five species of prairie dog found in North America and is listed as Threatened under the United States' Endangered Species Act and as Endangered by the IUCN Red List of Threatened Species (Roach, 2018). Due to heavy range-wide eradication campaigns during the 20th century, as well as ongoing habitat loss/fragmentation and epizootic outbreaks of plague, the Utah prairie dog has been reduced from a range-wide estimate of 95,000 individuals (estimated in the 1920s) to approximately 8969 today (total 2018 spring count of adults) (Collier & Spillett, 1973;Kavalunas & Day, 2018). In 1972, a recovery plan for the Utah prairie dog was enacted that focused on translocating individuals from private land to protected public lands (McDonald, 1993;United States Fish & Wildlife Service, 2012).
Today, Utah prairie dogs exist in small, isolated populations, making them less demographically stable and more susceptible to the erosion of genetic variation through genetic drift (Gilpin & Soulé, 1986;Wright, 1931). Further, Utah prairie dogs are highly social mammals that live in small family groups called coteries (consisting of one or two unrelated adult males, a group of related females, and their young; Hoogland, 2006), with adjacent coteries forming a colony.
In contrast to many mammalian species that exhibit natal dispersal, Utah prairie dogs rarely leave their natal coteries unless nearly all of the individuals in a coterie are gone (Hoogland, 2013), for example following an outbreak of plague. Any such dispersal that does occur is likely male-biased and to nearby coteries, often within the same colony (Hoogland, 2013). Due to prairie dog behavior coupled with the fragmented structure of their habitat patches, these colonies likely exhibit metapopulation dynamics (Brown et al., 2016).
Translocations are a common practice in the management of prairie dogs (United States Fish & Wildlife Service, 2012). This is an effective tool to combat the loss of genetic variation experienced through a history of population bottlenecks, eradication campaigns, disease, and limited natural gene flow. However, a large degree of variation in habitat exists across the Utah prairie dog species range.
Utah prairie dogs can be found from valley floors to montane habitats, with elevations ranging from roughly 1.5 to 2.9 km above sea level (Roach, 2018). Combined with a lack of natural gene flow, this habitat heterogeneity could encourage local adaptation (Blanquart et al., 2013). Under this scenario, translocating individuals from different areas of the species range could introduce maladaptive genes to other colonies, leading to outbreeding depression. Even though fears of outbreeding depression may be inflated Ralls et al., 2018), translocations between locally adapted populations could still have potentially disastrous consequences for the viability of prairie dogs. By characterizing not only neutral genetic markers, but also those under selection, we may avoid stifling the evolutionary potential of Utah prairie dogs.
The objective of this study was to use a population genomics approach to investigate the effect of gene flow, genetic drift, and divergent selection on the maintenance of genetic variation in Utah prairie dogs across their species range. To accomplish this objective, we carried out five aims. First, we characterized patterns of genetic structure and gene flow among Utah prairie dog populations. Second, we evaluated the impact of sex-biased dispersal on the maintenance of genetic variation by identifying differences in patterns of genetic structure and gene flow for females and males separately. Third, we characterized the impact of genetic drift on the erosion of genetic variation and genetic differentiation. Fourth, we identified loci under divergent selection and compared patterns of population divergence at these loci against neutral loci. Fifth, we examined how the environment might influence local adaptation by identifying genotype-environment associations (GEAs). Our genome-wide approach allows us to harness information in both neutral and adaptive loci to tailor conservation activities to maximize the success of recovery efforts without incurring the potentially substantial costs that could result from translocating locally adapted individuals.

| Generating the SNP dataset
We trapped Utah prairie dogs and collected hair and whiskers during a field trial of a sylvatic plague vaccine (Rocke et al., 2017).
Samples were collected in 2014 from two plots within each site (assumed to be a single prairie dog colony), with plots located in proximity (0.15-2.10 km). Due to a high degree of movement between plots observed during the vaccine field trial, plots were treated as a single site for our analyses. Individuals were sampled throughout the Utah prairie dog range at three sites near Cedar City and Panguitch, Utah (CCUT) and four high-elevation sites within the Awapa Plateau (HEUT; Figure 1). From the total estimated area of occupancy for Utah prairie dogs (2800 ha; Roach, 2018), we sampled 29.2 ha for CCUT sites and 54.7 ha for HEUT sites (Rocke et al., 2017). These sites are within the three recovery units ( (Catchen et al., 2011(Catchen et al., , 2013, following the proposed workflow outlined by Rochette and Catchen (2017 replicates. For STRUCTURE, we used the alternative prior for population-specific ancestry (α = 1) because we had unequal sampling among sites (Wang, 2017). In order to accurately resolve the number of genetic clusters (K) using STRUCTURE, we used a combination of the LnP(D) and Δ K as outlined in Janes et al. (2017) calculated using STRUCTURE HARVESTER v.0.6.93 (Earl & vonHoldt, 2012). We also adopted a 'hierarchical STRUCTURE analysis' approach where each genetic cluster was analyzed iteratively in a new STRUCTURE TA B L E 1 Measures of genetic variability in Utah prairie dogs (Cynomys parvidens) for (a) each of the sampled sites and (b) the two putative populations identified using both STRUCTURE and DAPC based on 2955 variable SNP loci run in order to gauge substructure (Vähä et al., 2007). We used the program CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007) to assign individuals to genetic clusters using q-values from STRUCTURE.
For the DAPC, the analysis was first performed unsupervised (no prior knowledge of groups) using the sequential K-means clustering algorithm executed through the find.clusters function in adegenet (Jombart & Ahmed, 2011;Jombart et al., 2010). We then performed the DAPC analysis supervised by using the Bayesian information criterion (BIC) to determine the final value of K.
To determine patterns of gene flow, we calculated the pairwise genetic differentiation between sites using Weir and Cockerham's F ST and Jost's D using the R packages hierfstat (Goudet & Jombart, 2015) and mmod (Winter, 2012), respectively. We tested for significant population divergence using both F ST and Jost's D with 1000 random permutations and corrected p-values for multiple comparisons using a false discovery rate (FDR) correction (used the 'p.adjust' function in R v3.5.2; Benjamini & Hochberg, 1995;R Core Team, 2018). To test for sex-biased dispersal, we repeated the above analyses on males and females separately. We estimated the effective number of migrants (N M ) based on the number of private alleles using the R package genepop (Barton and Slatkin, 1986;Rousset, 2008). We estimated relative migration rates among sampling sites using N M with the divMigrate function in the R package diveRsity and identified significant migration rates using 10,000 bootstrap iterations (Alcala et al., 2014;Keenan et al., 2013;Sundqvist et al., 2016). Migration networks were created using the R package qgraph (Epskamp et al., 2012). Values <1 indicate that populations will diverge due to drift.

| Characterizing the effect of genetic drift
Genetic drift erodes standing genetic variation where rare alleles face a greater chance of being lost due to random chance. As populations decrease in size and become isolated, the effects of genetic drift grow in significance, which accelerates the erosion of genetic variation. Thus, populations in which genetic drift is the primary driver of genetic structure are predicted to have (a) less genetic variation, (b) fewer private alleles (or number of unique alleles for a specific location), and (c) a higher degree of inbreeding than populations in which genetic structure is shaped by other mechanisms. We char- We estimated N e for each site using the linkage disequilibrium (LD) method implemented in the program NeEstimator v2.1 (Do et al., 2014). One assumption to calculate N e using the LD method is that linkage disequilibrium at independent loci in randomly mating, closed populations comes exclusively from genetic drift (Hill, 1981).
To meet this assumption, we created a neutral set of loci by excluding loci potentially under the influence of natural selection; those identified as outliers by either BayeScan or PCAdapt for this analysis (see next section for details).
We tested each of the seven prairie dog sampling sites as well as each genetic cluster for evidence of genetic bottlenecks using the program Bottleneck v.1.2.02 (Cornuet & Luikart, 1996;Piry et al., 1999). Like the estimation of N e , only neutral loci were used to test for evidence of bottlenecks. We used the infinite alleles model (IAM) and tested for significant heterozygosity excess compared to the level predicted under mutation-drift equilibrium using standardized differences tests (Cornuet & Luikart, 1996).

| Neutral loci versus loci under selection
Outlier loci were detected via a Bayesian method conducted in BayeScan 2.1 (Foll & Gaggiotti, 2008) and a nonconstrained ordination method executed in the R package PCAdapt (Luu et al., 2019). BayeScan is based on the multinomial-Dirichlet model and identifies differences in allele frequencies between subpopulations (we used sampling sites) and the common gene pool of all subpopulations (measured as a subpopulation specific F ST coefficient). For the BayeScan method of outlier detection, we defined populations as the sampling sites (n = 7) and used the default parameter settings (prior odds for neutrality = 10; odds are that 1 locus is under selection for every 10 neutral loci). We directly controlled the false discovery rate (FDR) by setting the target FDR (q-value) to limit the proportion of false positives to 5% (Foll & Gaggiotti, 2008). For the PCAdapt method, predefined populations are not required. We specified the K parameter (the number of PCs to retain) in PCAdapt as 5 based on where the amount of variation explained by each PC decreases (identified as the point of inflection in the PCAdapt 'scree plot'). We defined outlier loci as those with q-values <0.05, meaning that 5% or less of loci identified as outliers are potentially false positives. To generate the set of outlier loci, we retained loci that were identified as under divergent selection in both BayeScan and PCAdapt methods. We compared this outlier dataset to the neutral set of loci to assess the relative contribution of genetic drift and selection in shaping patterns of genetic structure. Specifically, we compared patterns of genetic structure for outlier and neutral loci using the program STRUCTURE and by using DAPC in adegenet (Jombart & Ahmed, 2011;Jombart et al., 2010). Gene flow was also estimated for each set of loci using pairwise Jost's D and F ST . To determine if geographic distance alone drives patterns of genetic structure, we tested for isolation by distance (IBD) in the full genetic dataset (neutral and outlier loci), neutral loci only, and outlier loci only using Mantel tests between genetic and geographic distance carried out in the R package adegenet (Jombart & Ahmed, 2011;Jombart et al., 2010). We calculated genetic distance using linearized pairwise F ST values (F ST /(1 − F ST )) and geographic distance was calculated as the pairwise Euclidean distance between sites (km) using the pointDistance function in the R package "raster" (Hijmans, 2019). Significance of Mantel tests was assessed based on 999 replicates.

| Selection and environmental associations
We conducted a GEA analysis using a multivariate ordination method, redundancy analysis (RDA), implemented in the R package 'vegan' Oksanen et al., 2019). This was was assumed if VIF >4). Using the RDA, we identified outlier loci and the environmental variables most associated with those loci (Capblancq et al., 2018;Forester et al., 2018). To prepare the genetic data for the RDA, we replaced missing data with the most common genotype across individuals and converted the full SNP dataset to allele counts. We used an analysis of variance (ANOVA) to quantify how well the full RDA model as well as each RDA axes explained genetic variation using the ANOVA.cca function in the R package vegan with 999 permutations (Legendre et al., 2011;Oksanen et al., 2019). Outlier loci were identified as those with loadings ± 3 standard deviations away from the mean loading. We used all of the significant constrained axes of the partial RDA for outlier detection (SNPs that load heavily on the axes are more likely to be under selection, p < 0.05). To infer how each SNP relates to the environmental variables, we identified the environmental variable that was most strongly correlated to each outlier SNP based on correlation coefficients. To illustrate these environment-SNP associations, we created triplots of all RDA axes using the R package vegan (Oksanen et al., 2019). For the triplots, symmetrical scaling (using the square root of eigenvalues) was used for both SNP and individual scores.

| Generating the SNP dataset
After initial filtering steps, a total of 3549 variable SNP loci were retained. An additional 594 loci were removed due to exceptionally high depth of coverage, indicating potential paralogs (2× the mode of the mean depth of coverage for each locus; mode = 17.55).
The final genomic dataset contained 2955 variable SNP loci with a mean depth of coverage of 20.08 for all individuals (mean per locus depth of coverage ranged from 7.68 to 35.58; Figure S3). Four individuals with a high amount of missing data (>30%) were removed from the dataset, leaving a total of 233 individuals used for subsequent analyses (13-58 individuals per site, mean = 33.29). These 233 individuals used for analyses had a mean of 4.2% missing data ( Figure S4).

| Characterizing patterns of genetic structure, gene flow, and sex-biased dispersal
Differentiation among sampling sites was observed using PCA ( Figure 3). The most supported number of genetic clusters (K) was two, using both the program STRUCTURE (based on Ln(K) and the ∆K method) and the unsupervised clustering method using DAPC ( Figures S5 and S6). In the two clusters, all individuals from the CCUT sites grouped to make one genetic cluster and all individuals from HEUT sites grouped to make the second genetic cluster.
When using a supervised clustering approach for DAPC, additional clustering solutions (K = 2-4) were informative for describing genetic structure based on BIC ( Figure S7). For example, under a K of four for DAPC, CCUT3 formed a separate genetic cluster and an additional genetic cluster included all individuals from HEUT2 and 11 individuals from HEUT4 ( Figure S8). This pattern was also observed using the hierarchical clustering approach in STRUCTURE ( Figure S6).
We found similar patterns of population differentiation with pairwise values of F ST and Jost's D (Table S2)

| Neutral loci versus loci under selection
Using BayeScan, we identified 754 loci as outliers, of which, 303 loci were under divergent selection. Using PCAdapt, 531 loci were identified as under selection ( Figure S9). For our outlier dataset, only loci that were identified as potentially under divergent selection in both BayeScan and PCAdapt were used (n loci = 51; 1.7% of total loci). For the neutral dataset, all loci as identified as an outlier in either program were removed (n loci = 1792; 60.6% of total loci).
We found significant patterns of isolation by distance using the full

| Selection and environmental associations
We found that a large proportion of the climatic variables were highly correlated (|r| > 0.70) to elevation. Therefore, we removed highly correlated variables as well as variables with low variation across sites. We reduced the four landscape (three landcover and elevation) and 19 climatic predictor variables to 4 total variables - Outlier (50.4%) followed by the second (27.5%), third (12.2%), and fourth (9.9%) axes. By using a ±3 standard deviation cutoff for SNP loadings on each RDA axis, we identified 141 outlier loci. Of the outlier loci, we found 46 SNPs were highly correlated with elevation, 2 with forests, 48 with precipitation seasonality, and 45 with temperature (Table S3). In the triplot, environmental variables (black arrows), SNPs (gray points), and individuals (colored circles) were arranged in ordination space based on their relationship with each ordination axis (Figure 7; Figure S11). From these plots, it was evident that temperature and precipitation seasonality loaded strongly on the first axis of the RDA and individuals from HEUT1 and CCUT3 showed strong associations with more precipitation and lower temperatures.
Conversely, individuals from HEUT2 and HEUT4 were associated with warmer temperatures during the driest quarter and lower precipitation seasonality. Individuals from HEUT4 were strongly associated with a higher proportion of forest within a 5 km buffer while HEUT1 was more associated with a lower proportion of forest. We then looked at SNP loadings with environmental variables to see how they relate in ordination space (Figure 8; Figure S12). Of the outlier SNPs, 10 were identified on the first RDA axis, 107 on the third axis, and 24 on the fourth axis ( Figure 8; Figure S12).

| D ISCUSS I ON
Through the use of high-resolution, genome-wide markers we have demonstrated that researchers can characterize patterns of genetic structure and gene flow, elucidate demographic processes, compare the effects of drift and selection on population divergence, and identify GEAs in a threatened species. This information is invaluable in tailoring the conservation policy of threatened species to mitigate the loss of genetic variation and maintain evolutionary potential. We found that Utah prairie dogs show high genetic structure, limited gene flow, and a lack of genetic variation across their range. This pattern can be explained sufficiently well by demographic processespopulation isolation and genetic drift. These processes seem to have had a greater impact on the CCUT sites, where genetic variation was lower than in the HEUT sites. However, selection also plays a role in the divergence of Utah prairie dog populations, primarily between low-elevation (CCUT) sites and high-elevation (HEUT) sites. to rarely leave natal colonies (Hoogland, 2013). This was reflected in both Jost's D and F ST metrics. Higher rates of migration between neighboring sites combined with strong evidence of isolation-by-distance indicate that when prairie dogs do disperse, it is likely to nearby colonies and that long-distance dispersal is rare or nonexistent.

Limited genetic variation in
These combined analyses lend support to Utah prairie dogs functioning as a metapopulation with the sites acting as subpopulations within the genetic clusters. Patterns of metapopulation dynamics, largely driven by plague epizootics, have also been documented in F I G U R E 7 Redundancy analysis (RDA) triplots showing environmental associations with outlier loci (RDA axes 1 and 2). The dark gray dots located at the center of the plot represent single nucleotide polymorphisms (SNPs), the colored points refer to individuals (colors represent which site they were sampled from), and black vectors represent environmental variables. SNP and individual RDA scores are scaled by the square root of their eigenvalues. The direction of the arrows indicates the correlation of the environmental variable with each axis black-tailed prairie dog (Cynomys ludovicianus) colonies (Antolin et al., 2006;Sackett et al., 2012). As plague also causes high mortality in Utah prairie dogs, it may serve to drive metapopulation dynamics in these populations as well. Further, additional colonies exist between CCUT and HEUT sites. These sites could act as "stepping stones" to facilitate dispersal that would result in a stronger pattern of IBD. We did not detect evidence of sex-biased dispersal, which we expected based on the observations in Hoogland (2013). This could be explained by a number of factors. First, young males are likely to remain in their natal coterie until they mature. If a large proportion of the males we sampled were juveniles, which was the case in another Utah prairie dog study (Brown et al., 2016), then our estimates of male dispersal could be unrealistically low. Male prairie dogs also actively defend their coterie (and females) and have even been found to occasionally engage in infanticide and cannibalism (Hoogland, 2007). This means it is possible that males disperse as observed by Hoogland (2013), but individuals that disperse experience high mortality by either other prairie dogs or by predators before they pass on their genes. Small N e was also reported in another study of Utah prairie dog using different markers and methodology (Brown et al., 2016). Drift acts more rapidly in small versus large populations (Lacy, 1987), which means that without gene flow we should expect genetic diversity to decline quickly in Utah prairie dog populations, despite positive demographic trends. Facilitated gene flow (i.e., translocation) from populations with unique alleles and high genetic variation would increase N e of small populations and slow drift-based erosion of diversity (Laikre et al., 2016). However, this strategy should be pursued with caution, as maladapted individuals introduced to new environments could reduce population fitness, as shown by the low establishment success of poorly adapted seeds (Kulpa & Leger, 2013) and reduced hatching success in non-native habitat in pike (Esox lucius) (Berggren et al., 2016).
Whereas restricted gene flow leaves genetic drift and inbreeding unmitigated, it also facilitates adaptation to local environments (Barrett & Schluter, 2008;Savolainen et al., 2013;Slatkin, 1987). The results from BayeScan and PCAdapt in this study reflect the observed high divergence and low rates of relative migration (N M ) between high-elevation (HEUT) sites and low-elevation (CCUT) sites.
The RDA analysis expanded on this pattern of divergence by showing that outlier loci were associated with high-elevation environments, a greater proportion of forests, colder temperatures during the driest quarter, and greater precipitation seasonality. Differences in habitat and climate between HEUT and CCUT sites indicate that local adaptation could occur, and our outlier and GEA analyses indicate that adaptation is ongoing. Other studies have also found that the HEUT sites harbor a greater diversity of flea and small mammal species compared to CCUT sites (Bron et al., 2018;Russell et al., 2018). This variation in environment and species composition likely plays a role in plague dynamics (Eads & Hoogland, 2017;Williams et al., 2013) that could affect both population viability as well as reduce the standing genetic variation required for adaptation to occur.
A large proportion of our SNP loci in our study were identified as under divergent selection (10% using BayeScan and 18% using PCAdapt). This is somewhat larger than the proportion of SNP loci identified in studies of other mammals (7% of loci using PCAdapt in zebras (Equus quagga; Pedersen et al., 2018); 4.4% loci for BayeScan and 6.7% of loci for PCAdapt in red foxes (Vulpes vulpes; Roberts, 2019)), and likely reflects the elevated opportunity for adaptation in a species with high site fidelity, small populations, and short generation time. Genome scans can be prone to false positives and be challenging to interpret when population structure is present (Capblancq et al., 2018;Forester et al., 2016Forester et al., , 2018Hoban et al., 2016). We aimed to minimize these limitations by generating a conservative set of criteria for outlier detection (those that were identified using two different approaches; 1.7% of total loci) and reducing the potential for false positives through FDR control. Further, comparing environmental variables to allele counts in a constrained ordination framework (RDA), as was done in this study, has low false-positive rates and high true-positive rates across demographic history and sampling schemes (Capblancq et al., 2018;Forester et al., 2016Forester et al., , 2018. Denser sampling across the genome, for example incorporating whole-genome sequences, would provide a more accurate representation of genomic variation across all linkage blocks in order to better gauge local adaptation and identify loci under selection (Hoban et al., 2016;Lowry et al., 2017). Further, experimental work such as common garden or reciprocal transplant experiments would help to validate true connections between the genotypes and their associated environmental variables.

| CON CLUS IONS
With continued change in land use and climate, maintaining genetic variation is of paramount importance for the survival of threat- Translocations initiate gene flow and increase the genetic variation in the recipient population (Whiteley et al., 2015), and have been successful in threatened species conservation (Bateson et al., 2014;Johnson et al., 2010). However, translocations could exacerbate population declines in threatened species by introducing maladaptive genes into populations that are locally adapted (outbreeding depression; Frankham, 2015;Weeks et al., 2011).
Incorporating selection directly into translocation strategies could help to mitigate the risk of translocation-mediated breakdown of local adaptation while increasing standing variation (Flanagan et al., 2018;Harrisson et al., 2017). We also suggest that Utah prairie dogs are low risk of outbreeding depression based on two of the three screening criteria outlined by Frankham et al. (2011) populations are chromosomally compatible, and isolation of populations occurred within the last 500 years. The third criteria for low outbreeding depression risk, that populations are not adapted to vastly different habitats, is indirectly supported by our data indicating natural gene flow between high and low-elevation populations (HEUT sites and CCUT3) and an overall pattern of IBD.
With low levels of genetic variation, especially in the CCUT sites, and limited gene flow species-wide, the consequences of genetic drift in Utah prairie dogs may outweigh concerns of outbreeding depression from moving potentially locally adapted individuals (Hereford, 2009;Ralls et al., 2018;Weeks et al., 2016). Specifically, translocations from higher-diversity HEUT sites to CCUT sites may serve to increase genetic variation in CCUT sites. The effect of translocations on local adaptation in CCUT is not easily predicted, but recent evidence (Fitzpatrick et al., 2020) indicates that gene flow into small, isolated populations can bolster overall genetic variation without erasing local adaptation, an ideal management outcome for Utah prairie dogs.

ACK N OWLED G M ENTS
We thank Loren Cassin Sackett for providing the Gunnison's prairie dog genome that was used for the assembly of ddRAD sequences.
We thank Texas A&M for performing all the sequencing used in this study. We thank Carly Malave for DNA extraction from all the hair and whiskers used in this study. We thank Adam Kavalunas and Robin Russell for reviewing this manuscript and providing helpful comments. We thank the Associate Editor and the anonymous reviewers for their suggestions for revisions that improved the manuscript. Funding was provided by the Morris Animal Foundation (D14ZO-405 and D14ZO-044). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genomic data used in this study are available from the Dryad Digital