The genetic legacy of 50 years of desert bighorn sheep translocations

Abstract Conservation biologists have increasingly used translocations to mitigate population declines and restore locally extirpated populations. Genetic data can guide the selection of source populations for translocations and help evaluate restoration success. Bighorn sheep (Ovis canadensis) are a managed big game species that suffered widespread population extirpations across western North America throughout the early 1900s. Subsequent translocation programs have successfully re‐established many formally extirpated bighorn herds, but most of these programs pre‐date genetically informed management practices. The state of Nevada presents a particularly well‐documented case of decline followed by restoration of extirpated herds. Desert bighorn sheep (O. c. nelsoni) populations declined to less than 3,000 individuals restricted to remnant herds in the Mojave Desert and a few locations in the Great Basin Desert. Beginning in 1968, the Nevada Department of Wildlife translocated ~2,000 individuals from remnant populations to restore previously extirpated areas, possibly establishing herds with mixed ancestries. Here, we examined genetic diversity and structure among remnant herds and the genetic consequences of translocation from these herds using a genotyping‐by‐sequencing approach to genotype 17,095 loci in 303 desert bighorn sheep. We found a signal of population genetic structure among remnant Mojave Desert populations, even across geographically proximate mountain ranges. Further, we found evidence of a genetically distinct, potential relict herd from a previously hypothesized Great Basin lineage of desert bighorn sheep. The genetic structure of source herds was clearly reflected in translocated populations. In most cases, herds retained genetic evidence of multiple translocation events and subsequent admixture when founded from multiple remnant source herds. Our results add to a growing literature on how population genomic data can be used to guide and monitor restoration programs.

examples where translocations have restored populations (reviewed by Seddon & Armstrong, 2016). The successful reintroduction of gray wolves (Canis lupus L.) into the Greater Yellowstone Ecosystem led to the re-establishment of ecosystem processes (Ripple & Beschta, 2012), while the creation of offshore populations of New Zealand's South Island saddleback (Philesturnus carunculatus Gmelin, 1789) prevented near-certain extinction in response to invasive rodents (Hooson & Jamieson, 2003). These  Translocation activities can have considerable impacts on the genetics of natural populations, including reduction in genetic variation, erosion of local adaptation, and changes to preexisting landscape genetic structure (Laikre, Schwartz, Waples, & Ryman, 2010).
However, translocations involving organisms with strong site fidelity pose unique difficulties, as these populations are naturally prone to the erosion of genetic diversity when physically isolated from other populations (Segelbacher, Höglund, & Storch, 2003;Westemeier et al., 1998).
The most widely cited estimate of the historical number of bighorn sheep in western North America before extensive European settlement is ~1.5-2 million individuals (Seton, 1929). Dramatic declines occurred throughout the late 1800s and early 1900s, largely due to unregulated hunting, overgrazing, and susceptibility to diseases often transmitted from domestic sheep (Ovis aries L.) and goats (Capra hircus L.) (Buechner, 1960;Cassirer et al., 2018;Valdez & Krausman, 1999). Declines ultimately resulted in the extinction of at least one taxonomically contentious evolutionary lineage (Badlands bighorn sheep;O. c. auduboni Merriam, 1901) and perhaps other unrecognized lineages (Malaney et al., 2015). In the United States, less than 20,000 individuals remained by 1960 (Buechner, 1960), with most individuals found in small, isolated herds scattered throughout the remaining portion of the range. Several genetic consequences likely stemmed from these declines, including increased isolation of fragmented populations and reduced genetic diversity within herds.
In an effort to repopulate previously occupied mountain ranges, augment genetic variation within isolated remnant herds, and increase connectivity among populations, state and federal agencies developed management programs that conducted several hundred translocations throughout the late 1900s (Wild Sheep Working Group, 2015). As part of this enterprise, the Nevada Department of Wildlife (NDOW) undertook an extensive series of translocations that have spanned five decades (1968-present), raising the statewide population estimate from less than 3,000 individuals in 1960 to at least 12,000 individuals in 2018 (NDOW, 2001, NDOW unpublished data). NDOW's translocation strategy was guided by Cowan (1940) and Hall's (1946) Table S1; Malaney et al., 2015;Wakeling, 2015;Wild Sheep Working Group, 2015). It is important to note that Wehausen andRamey (1993, 2000) subsequently disputed Cowan (1940) and Hall's (1946) hypothesis and instead suggested that desert bighorn sheep historically occupied all of Nevada, with a small-horned Great Basin lineage in the north and a large-horned Mojave lineage in the south (see Figure 1). Thus, there is uncertainty about the extent to which translocated herds may have been locally adapted to different environmental conditions, as well as the consequences of translocation history for population genetic diversity and structure across the landscape.
Here, we used genotyping-by-sequencing (GBS) to generate population genomic data for desert bighorn sheep sampled from the Great Basin and Mojave Deserts to explore the population genetic context and consequences of these translocations. Our goals were to: (a) characterize genetic diversity and differentiation among herds, especially among Mojave source herds; (b) evaluate patterns of admixture in reintroduced herds in light of translocation histories; and (c) quantify the degree of genetic diversity within translocated herds composed of either single or mixed origins relative to source populations. Our results provide a reference for continued translocation decisions and a baseline for understanding how past and future population responses might relate to genetic variation within and among herds. More broadly, our study illustrates how high-throughput sequencing approaches can be used to illustrate the population genetic context and consequences of translocation activities spanning several decades.

| Sample DNA collection
We obtained muscle, liver, or blood samples from a hunter harvest return program and a long-term herd-monitoring program (NDOW).
We maximized spatial coverage of samples to represent as many  Table S1 for the full translocation history of desert bighorn sheep in Nevada. Black circles represent the location of the two major cities in Nevada: Reno in the northwest and Las Vegas in the south. The geographic bounds of the Mojave and Great Basin Deserts (purple and green, respectively) were taken from the United States Environmental Protection Agency level III ecoregions (Omernik & Griffith, 2014) extant herds as possible, but actively monitored populations or those with higher harvest rates were disproportionally represented.
DNA was extracted using Qiagen DNeasy Blood and Tissue Kits and quantified via spectrophotometry with a QIAxpert device (Qiagen Inc., Valencia, CA).

| Genotyping-by-sequencing
Reduced-representation libraries for Illumina sequencing were prepared using a GBS protocol (Parchman et al., 2012) analogous to ddRADseq (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012). First, extracted DNA was digested using two restriction enzymes with six base pair (bp) recognition sites (EcoRI and MseI). Raw sequencing data were filtered for contaminant DNA (e.g., PhiX, E. coli), low-quality reads, and Illumina adaptors using bow-tie2 _ db (Langmead & Salzberg, 2012) and a suite of bash and Perl scripts designed for this purpose (https://github.com/ncgr/tapioca). Variant sites (i.e., single nucleotide polymorphisms; SNPs) were identified, and genotype likelihoods were calculated using SAMtools v1.3 and BCFtools v1.3 . By using genotype likelihoods rather than calling genotypes categorically, these estimates reflect genotype uncertainty arising from variation in sequencing depth across individuals and loci (e.g., Nielsen, Paul, Albrechtsen, & Song, 2011). We set a minimum base quality of 20, a minimum mapping quality of 20, and required variants to have a minimum site quality of 20 and minimum genotype quality of 10. Final quality filtering of variant sites was performed with VCFtools v0.1.14 (Danecek et al., 2011), and only biallelic SNPs were retained. Loci were only included in the final dataset if their minor allele frequency (maf) was greater than 0.05 and if at least 60% of individuals had at least one read present at the locus. Additionally, individual sheep were removed if they had >50% missing data. Finally, all loci with a mean individual coverage greater than 10X were removed in an effort to filter loci that could represent mis-assembly of paralogous regions. Considering our sequencing data represent low to medium levels of coverage, we used probabilistic methods to infer genotype probabilities while accounting for stochastic variation among loci and individuals in coverage and quality (see below).  (McQuivey, 1976b(McQuivey, , 1976c, which we hereafter refer to as the Muddy Mountains population for simplicity. Each of the remnant source populations is thought to have been demographically stable throughout the mid-1900s (McQuivey & Leslie, 1976McQuivey, 1976aMcQuivey, , 1976cMcQuivey, , 1979, though the Mormon Mountains suffered a single, rapid die-off (>50% individual mortality) following a population expansion in the early 1980s (McQuivey, 1982).

| Population genetic analyses
Individuals from source herds have been used to re-establish more than half of the desert bighorn sheep herds in Nevada, as well as founding other populations in Colorado, Texas, and Utah (Wakeling, 2015;Wild Sheep Working Group, 2015).
We used a hierarchical Bayesian model (entropy;  that is based on the correlated allele frequency model of structure (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) to infer the number of ancestral genetic clusters (k), estimate ancestry coefficients for each individual, and estimate genotype probabilities for each individual at each locus. Importantly, this model utilizes a population allele frequency prior and incorporates information about genotype uncertainty arising from variation in sequencing coverage, sequence error, and alignment error during parameter estimation .
In order to determine the most probable k in our dataset, we conducted five replicate analyses of entropy from k = 2 to k = 7 for the full dataset and k = 2 to k = 4 for the remnant source dataset. The fit of each model was evaluated using a deviance information criterion (DIC), with smaller DIC values consistent with better model fit. Analyses for k > 7 in the full dataset and k > 4 in the source dataset failed to consistently converge and had high DIC values, so they were not considered further. To provide MCMC sampling with a starting point and to facilitate chain convergence, we conducted principal component analysis (PCA) on a covariance matrix of the genotype likelihoods and used k-means clustering based on principal components (PCs) to cluster individuals. Starting membership probabilities to the clusters were then calculated for each individual with linear discriminant analyses (LDA) on all PCs. Both k-means clustering and LDA were performed using the MASS package (Venables & Ripley, 2002) Team, 2017). entropy analyses were based on 100,000 MCMC iterations with a burn-in of 30,000 and thinning every tenth step. We additionally used PCA to summarize genetic variation among individuals and identify genetic structure among herds. The genotype covariance matrix from genotype probabilities generated with entropy was used as the input for PCA, which was performed using the prcomp function in R. The matrix of genotype probabilities is available at Dryad (https://doi.org/10.5061/dryad.25f502n).
To quantify differentiation among hunt units, Hudson's F ST (Hudson, Slatkin, & Maddison, 1992) Table 1). Haversine geographic distances were calculated using the fossil package (Vavrek, 2011) in R based on the midpoints of latitude and longitude for each hunt unit (Table S2).

| Genetic structure among the remnant, source herds
We used a hierarchical Bayesian model (entropy) to investigate potential fine-scale genetic structure among four source herds (units 212, 268, 269, 271) that have never received translocated individuals (Table   S1). Based on DIC, the model with four genetic clusters (k) was better supported than models with k = 2 or k = 3 (Table S3) Table 2).

| Genetic consequences of translocations
Comparison of the entropy models using the entire dataset of 303 individuals suggested that models from k = 2 to k = 6 had roughly equivalent support (Table S4). Qualitatively, these models were complementary to one another, consistent with the presence of hierarchical genetic structure. k = 6 had the lowest mean DIC and k = 4 had the highest mean DIC, but the magnitude of difference between all models was fairly low (Table S4)

| Fine-scale genetic structure among remnant Mojave populations
Our results illustrate fine-scale population genetic structure among remnant populations of desert bighorn sheep in the northern TA B L E 2 Pairwise F ST (Hudson et al., 1992) (Figure 3; see Supplemental Results). Our results suggest evidence of population genetic structure across finer spatial scales than most past genetic studies of North American wild sheep, which detected genetic structure across various spatial scales using a broad range of molecular markers (e.g., Ramey, 1995;Luikart & Allendorf, 1996;Fitzsimmons, Buskirk, & Smith, 1997;Boyce, Ramey, Rodwell, Rubin, & Singer, 1999;Gutiérrez-Espeleta, Kalinowski, Boyce, & Hedrick, 2000;Worley et al., 2004;Miller, Poissant, Kijas, & Coltman, 2011;Buchalski et al., 2015Buchalski et al., , 2016Kardos et al., 2015;Malaney et al., 2015;Sim, Hall, Jex, Hegel, & Coltman, 2016). This increased spatial resolution could be a product of the geographic distribution of isolated mountain ranges that bighorn sheep occupy in Nevada (Figure 1), but was also likely influenced by the relatively large number of markers we employed (~17,000 SNPs) compared to most past studies.
The relatively fine geographic scale at which we detected genetic structure is likely associated with both the natural history of bighorn sheep and the recent history of human activity in this region. First, the life history of bighorn sheep predisposes populations to genetic differentiation, as individuals often exhibit high site fidelity to natal habitats with access to high-quality forage, escape terrain, and water sources (McQuivey, 1978). Furthermore, the skewed mating ratio of bighorn sheep, where a few rams account for most of the successful mating events, (Coltman et al., 2002;Hogg & Forbes, 1997;Martin et al., 2016;Pelletier et al., 2006) can lower effective population sizes, intensify genetic drift within herds, and lead to population structure across increasingly fragmented landscapes, as seen in other organisms with polygynous mating systems (e.g., Coltman, Pilkington, & Pemberton, 2003;Bouzat & Johnson, 2004;Shafer, Côté, & Coltman, 2011;Jahner et al., 2016;Dotsev et al., 2018). isolating the River Mountains herd from other populations (Leslie, 1977). Subsequently, a number of artificial water sources (i.e., guzzlers) were installed in the River Mountains (Leslie & Douglas, 1979;McQuivey & Leslie, 1976), minimizing the need for individuals to travel for water. Despite past reports suggesting that individuals migrated between the Muddy and River Mountains (Denniston, 1965), our results suggest substantial gene flow has not recently occurred between these herds ( Figure 2). Thus, these naturally differentiated herds could have been further isolated and subdivided by the past 80 years (~11 sheep generations) of human activity and infrastructure. It is worth noting, however, that the effectiveness of anthropogenic barriers to sheep dispersal can vary substantially over time (Epps et al., 2018), and future changes could influence the degree of differentiation between these two remnant populations.

| Remnant Great Basin genetic ancestry
The original strategy for repopulating bighorn sheep to the Great Basin Desert was guided by a contentious taxonomic hypothesis (Cowan, 1940;Hall, 1946). Our results could lend support to Wehausen andRamey's (1993, 2000 ) alternative hypothesis of a historically widespread desert bighorn sheep range, with cranial morphology and population substructure strongly matching ecotypic differences among western North American deserts (i.e., a small-horned Great Basin lineage and a large-horned Mojave lineage). The herd we sampled on Lone Mountain (unit 212) has one of F I G U R E 3 Genetic structure of remnant and reintroduced herds of desert bighorn sheep (Ovis canadensis nelsoni) based on 17,095 SNPs. (a) For each individual, an ancestry coefficient was estimated for each of four genetic ancestries (k) with entropy   Table S1 for a full translocation history of all hunt units) the smallest mean horn sizes in Nevada (M. Cox, personal observation) and appears to be resistant to pneumonia despite the recent detection of Mycoplasma ovipneumoniae in the population (NDOW, 2017). This is also the most genetically differentiated population we sampled (Table 2, Figure S6) and may thus be a relict of an unrecognized Great Basin lineage of desert bighorn sheep. Nonetheless, the magnitude of differentiation between Lone Mountain and the Mojave source herds is modest (Table 2), so it is possible that the remnant Great Basin populations are simply recently differentiated herds that may be locally adapted to an ecosystem that differs dramatically from the Mojave Desert in both climate and vegetation (Beatley, 1975;Pavlik, 2008).  (Malaney et al., 2015;Wehausen, 1991).
For example, North Dakotan bighorn sheep had greater recruitment and projected population growth rates when they were sourced from environmentally similar populations (Bleich, Sargeant, & Wiedmann, 2018), lending support to recent calls for managers to consider using ecologically similar source populations for translocations (e.g., Lawrence & Kaye, 2011;Malaney et al., 2015;Biebach et al., 2016;Kronenberger et al., 2017). Although inferences regarding local adap-  Table 1). There was marginal evidence for elevated diversity in multiple source herds relative to single source herds ( Figure S8), which has also been reported for California bighorn sheep populations in Oregon (Olson, Whittaker, & Rhodes, 2013  . While this strategy is less than ideal for managing fragmented populations that are at elevated risk of losing genetic diversity (Frankham et al., 2017), the shortterm risks of disease transmission are currently too high to justify attempts to mitigate the more long-term risks of reduced genetic diversity. Looking forward, future success in this system will depend on how well translocations are leveraged to maintain genetic variation and maximize population persistence while preserving the identities of multiple evolutionary ancestries.

| Informing restoration with population genomic variation
Even prior to the advent of modern sequencing technologies, a pivotal question in restoration genetics was whether translocation pro- Thus, even though genetic rescue has been demonstrated to be an effective tool for restoring genetic diversity in small populations of organisms (Frankham, 2015), including bighorn sheep (e.g., Hogg, Forbes, Steele, & Luikart, 2006;Gompert, 2012;Miller, Poissant, Hogg, & Coltman, 2012;Olson, Whittaker, & Rhodes, 2012), more widespread genetic sampling of potential source populations prior to translocation may be an effective approach.
Although several thousand translocations have been conducted over the past century without genetic monitoring (Laikre et al., 2010), much insight could be gained by analyzing the population genetic context and consequences of restoration programs. Such analyses can provide a survey of genetic variation across candidate source populations, illustrate how translocation activities alter landscape genetic variation, and provide insight into evolutionary history that may be relevant for understanding local adaptation. In general, large-scale translocations will greatly impact the landscape genetic structure of populations relative to preexisting natural conditions. For example, a study of 72 lake trout (Salvelinus namaycush Walbaum 1792) populations found that natural lakes had lower genetic diversity and higher genetic differentiation than stocked lakes (Valiquette, Perrier, Thibault, & Bernatchez, 2014). Similarly, in our study, reintroduced populations lacked any signature of isolation by distance (Fig. S7b) and pairwise F ST values among remnant source populations were elevated relative to those among re-established populations (Fig. S6). Higher density population genetic data can also reveal previously unrecognized genetic structure (e.g., the Lone Mountain herd in our study) that could represent unique patterns of evolutionary history potentially associated with local adaptation. The preservation of locally adapted populations is a primary goal of many restoration efforts (McKay et al., 2005;Weeks et al., 2011), and while the presence of population genetic structure alone does not provide evidence for local adaptation, managers often consider the preservation of genetically differentiated populations with unique phenotypic variation. Finally, genetic data can identify natural corridors that allow for individual movements among populations that can be important for maintaining genetic diversity without human assistance (e.g., the pioneering herds in Figure 1; Epps, Wehausen, Palsbøll, & McCullough, 2010;Gilbert-Norton, Wilson, Stevens, & Beard, 2010). As the generation of high-throughput DNA sequencing data has become rapid and cost-effective, the analysis of such data could increasingly be used to guide and assess the restoration of populations.

| Conclusions
Emerging technologies that allow for increases in the extent of genomic sampling are now reshaping our understanding of the genetic consequences of conservation actions, even for organisms with complicated management histories (Allendorf, Hohenlohe, & Luikart, 2010;Shafer et al., 2015). However, relatively few studies have utilized such datasets to investigate the genetic legacy of reintroduction programs (e.g., Campbell, Kamphaus, Murdoch, & Narum, 2017;Grossen et al., 2018). By employing a GBS approach to investigate the distribution of desert bighorn sheep genetic variation across the Great Basin and northern Mojave Deserts, we uncovered an intricate genetic landscape structured by 50 years of translocation decisions. Furthermore, our results support the possible existence of a previously hypothesized Great Basin lineage of desert bighorn sheep that may require revised management consideration.
The population genetic patterns identified here will serve as a baseline that can be used to inform future translocation decisions, as well as a reference to understand how future population responses to disease outbreaks, climate change, and other environmental challenges are affected by genetic diversity and variation.