Genotyping‐in‐Thousands by sequencing reveals marked population structure in Western Rattlesnakes to inform conservation status

Abstract Delineation of units below the species level is critical for prioritizing conservation actions for species at‐risk. Genetic studies play an important role in characterizing patterns of population connectivity and diversity to inform the designation of conservation units, especially for populations that are geographically isolated. The northernmost range margin of Western Rattlesnakes (Crotalus oreganus) occurs in British Columbia, Canada, where it is federally classified as threatened and restricted to five geographic regions. In these areas, Western Rattlesnakes hibernate (den) communally, raising questions about connectivity within and between den complexes. At present, Western Rattlesnake conservation efforts are hindered by a complete lack of information on genetic structure and degree of isolation at multiple scales, from the den to the regional level. To fill this knowledge gap, we used Genotyping‐in‐Thousands by sequencing (GT‐seq) to genotype an optimized panel of 362 single nucleotide polymorphisms (SNPs) from individual samples (n = 461) collected across the snake's distribution in western Canada and neighboring Washington (USA). Hierarchical STRUCTURE analyses found evidence for population structure within and among the five geographic regions in BC, as well as in Washington. Within these regions, 11 genetically distinct complexes of dens were identified, with some regions having multiple complexes. No significant pattern of isolation‐by‐distance and generally low levels of migration were detected among den complexes across regions. Additionally, snakes within dens generally were more related than those among den complexes within a region, indicating limited movement. Overall, our results suggest that the single, recognized designatable unit for Western Rattlesnakes in Canada should be re‐assessed to proactively focus conservation efforts on preserving total genetic variation detected range‐wide. More broadly, our study demonstrates a novel application of GT‐seq for investigating patterns of diversity in wild populations at multiple scales to better inform conservation management.

Anthropogenic disturbances such as residential and commercial development, roads, agriculture, and forestry can lead to the isolation of large areas of formerly intact habitat that becomes disconnected over time. As a result, populations often are reduced in size and exhibit low connectivity and restricted gene flow (Pease, Lande, & Bull, 1989). Without sufficient genetic variation, populations of a species may be limited in their ability to adapt in the face of stochasticity, making them more vulnerable to extinction (Frankham, 2003). Over time, the extirpation of multiple isolated subpopulations can contribute to overall species decline in population size and erosion of genetic diversity. As a result, species existing in isolated populations often are of elevated conservation concern (Ehrlich, 1988).
Population genetic studies play an important role in understanding connectivity among geographically isolated populations, especially for rare or at-risk species. In particular, genetic assessments can help delineate conservation units, which are fundamental for enforcing legislation and prioritizing conservation resources (Funk, Mckay, Hohenlohe, & Allendorf, 2012;Green, 2005).

Multiple concepts for designating conservation units have been
proposed within various contexts, including evolutionarily significant units (ESUs) based upon historical patterns of genetic divergence between populations (Crandall, Bininda-Emonds, Mace, & Wayne, 2000;Ryder, 1986;Waples, 1991) and management units (MUs) that consider contemporary patterns of genetic differentiation among populations (Moritz, 1994). In Canada, designatable units (DUs) have been adopted by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) to define biologically relevant units for conservation below the species level (Green, 2005). DUs are determined based upon population discreteness and evolutionary significance, two criteria that can be examined using genetic data. According to COSEWIC, population(s) may be considered discrete if there is evidence for genetic distinctiveness, natural geographic disjunction, or occurrence in geographically varied ecological zones. Once populations are deemed discrete, they can then be evaluated for significance to determine whether they are integral to the evolutionary trajectory of the species, where unique genetic diversity would not be naturally restored if lost (Environment & Climate Change Canada, 2017). Measures of neutral genetic variation and gene flow among populations can provide evidence for discreteness, while the evolutionary significance of populations can be informed by estimating patterns of adaptive divergence (Allendorf, Hohenlohe, & Luikart, 2010;Funk et al., 2012).
One species of concern in Canada that requires an explicit evaluation of designatable unit status is the Western Rattlesnake (Crotalus oreganus). Although broadly distributed throughout the western half of the United States, this snake is classified as "threatened" by the Canadian Species At Risk Act and protected under the British Columbia (BC) Wildlife Act (COSEWIC, 2015). The area of occupancy of this species in Canada is restricted to BC and estimated at <1,500 km 2 (COSEWIC, 2015; Figure 1), representing one of the northernmost extensions for any viper in North America (Family Viperidae). At these northern latitudes, Western Rattlesnakes use traditional denning sites that enable overwinter survival (Charland, Nelson, & Gregory, 1993). The numbers of animals using these dens range from single digits to several hundred individuals, including adults and juveniles of both sexes (Winton, Taylor, Bishop, & Larsen, 2018). Following emergence from the dens in the spring, adult snakes migrate various distances to forage and mate before returning to the same dens in the fall (Gomez, Larsen, & Gregory, 2015; Harvey & Larsen, in press). Gravid females do not travel as far, remaining closer or even at den sites to give birth to their litters in late summer (Eye and Larsen unpubl.). Presently, the entire distribution of the Western Rattlesnake in Canada is managed as a single DU (COSEWIC, 2015), although five putative geographic regions have been identified (Figure 1). These regions are reasonably distinct, given that suitable habitat at this latitude is relatively scarce and naturally interrupted by geographic barriers. Moreover, anthropogenic modifications to the landscape, such as habitat loss and fragmentation by roads, likely constitute additional barriers to gene flow in this species (Maida et al., 2018;Winton et al., 2018).
In addition to the influence of natural and anthropogenic barriers to gene flow, range position can impact the extent and distribution of population genetic variation in broadly distributed species.
The central-marginal hypothesis (Eckert, Samis, & Lougheed, 2008) predicts that peripheral populations, particularly those in discontinuous patches, are expected to have lower genetic diversity within and greater differentiation among populations than those found toward the core of a species range. Despite smaller population sizes and an increased risk of extirpation, these peripheral populations may also have high evolutionary potential given they are likely subject to more heterogeneous environments than found in the range core (Cassel-Lundhagen et al., 2009;Lesica & Allendorf, 1995). Such conditions may lead to the development of unique traits and divergent genetic variation within peripheral populations that permit adaptation to environmental change (Lesica & Allendorf, 1995). For Western Rattlesnakes, a fuller understanding of the distribution of genetic variation within and among populations at the northernmost extent of their species range is warranted for evaluating discreteness, inferring significance, and refining conservation strategies.
In the current study, we used genotyping-by-sequencing (GBS) to characterize the distribution of genetic variation within and among Western Rattlesnake populations in Canada and neighboring Washington (USA). Given the protected status and venomous nature K E Y W O R D S designatable units, genotyping-by-sequencing, GT-seq, hibernacula, reptile, snake of the animal, we largely employed minimally invasive sampling (MIS) using cloacal swabs, as it reduced handling time and limited disturbance to individual snakes. Although MIS is often suboptimal for use with many GBS approaches (e.g., RADseq; Andrews, Barba, Russello, & Waits, 2018;Russello, Waterhouse, Etter, & Johnson, 2015), we recently demonstrated the utility of a multiplex amplicon sequencing approach, Genotyping-in-Thousands by sequencing (GT-seq; Campbell, Harmon, & Narum, 2015), for cloacal swab samples and other MIS sources (e.g., skin sheds) for the Western Rattlesnake (Schmidt, Campbell, Govindarajulu, Larsen, & Russello, 2020). Here, we used genotypic data collected from 591 individuals to infer patterns of hierarchical genetic structure within and among the five presumed geographic regions of Western Rattlesnakes in British Columbia (Canada) as well as in Washington (USA). We also examined fine-scale patterns of genetic diversity and connectivity among complexes of dens within regions. Lastly, we used the evidence collected to explicitly address DUs for Western Rattlesnakes in Canada and discuss the implications of these results for informing prioritization and management strategies.

| Study area and sample collection
Between 2015 and 2017, a total of 591 individuals were sampled, including from blood (n = 30), cloacal swabs (n = 552), buccal swabs F I G U R E 1 The distribution of den complex locations of Western Rattlesnakes across five main geographic regions in Canada, including Thompson-Nicola (TNI), Vernon (VER), Okanagan-Similkameen (OKS), Midway (MW), Grand Forks (GF), as well as Washington (WA), USA. Gray shading represents the current known distribution in Canada and Washington. Black-dotted circles encompass geographic regions defined by COSEWIC (2015); WA was also considered a separate geographic region. Colored pie charts indicate results of the hierarchical Bayesian clustering analyses for individuals within the 11 den complexes across the six geographic regions: TNI-W (n = 57), TNI-C (n = 35), TNI-E (n = 17), VER-CB (n = 30), VER-KL (n = 57), WL (n = 99), OS (n = 55), MW (n = 23), GF (n = 31), WA-N (n = 39), WA-S (n = 18). Each color within a pie chart represents the proportion of identity to each unique genetic cluster. The black arrow is pointing to the City of Kelowna, British Columbia, Canada (gray star). Gray shading in the figure inset displays the entire distribution of Western Rattlesnakes in North America with the approximate distribution in this study outlined in dark gray [spatial data from NatureServe and the International Union for All DNA samples obtained from live snakes were taken with the anterior portion of the animals restrained in a plexiglass handling tube (Murphy, 1971). Blood samples were collected by inserting a sterile syringe into the caudal vein of a snake, drawing 0.1 ml of blood, and then directly transferring the sample to 1 ml of Longmire lysis buffer (Longmire, Maltbie, & Baker, 1997) for preservation at 4°C. Cloacal samples were collected by gently inserting a sterile cotton-tipped swab into the cloacal vent, located on the ventral side of the snake near the tail, and then slowly rotating it 5-10 times. For juvenile and neonate snakes that were too small for cloacal sampling, we employed buccal swabs. To collect buccal cells, a cotton swab was gently inserted into the mouth of a snake and slowly rotated for 5-10 s, with care taken to avoid getting the swab stuck on the fangs. As with blood samples, buccal and cloacal swabs were then preserved in 1 ml of Longmire lysis buffer (Longmire et al., 1997) at 4°C. Shed skins were stored at −80°C as recommended by Burbrink and Castoe (2009).
In addition to collecting DNA samples, we also recorded the snout-to-vent length (SVL), sex, weight (g), and age class for individual snakes when possible. For the present study, we defined adults as those snakes with a SVL measurement of >540 mm either at the time of sampling or postsampling, based upon recent estimates for sexually mature Western Rattlesnakes (Maida et al., 2018). All

| DNA extraction
We extracted whole genomic DNA from the blood samples using a DNeasy Blood and Tissue Kit (Qiagen) following the manufacturer's protocols (including RNAse A treatment). Only slight variations were applied for the shed, buccal, and cloacal swab samples. For the shed skins, ~160-180 mg of each shed underwent a surface clean with a 3% bleach solution. After bleach treatment, we then washed the sheds with ethanol and left them to air-dry completely. Once dry, we manually pulverized each sample and placed them into separate 2-ml tubes. Proteinase K and lysis buffer were added to each tube, and samples were incubated for 48 hr at 56°C. The rest of the extraction procedure followed the manufacturer's protocol (including RNAse A treatment).
For the cloacal and buccal swab samples, we extracted DNA from the lysis buffer in which the swabs were stored. 180 μl of buffer from each sample was incubated with 20 μl of Proteinase K for 1 hr at 56°C. The rest of the extractions then followed the manufacturer's protocol (including RNAse A treatment). Once eluted, we quantified all DNA extracts using a Qubit™ 3.0 Fluorometer and Qubit™ dsDNA High Sensitivity Assay kit (Invitrogen).

| GT-Seq library preparation and genotyping
We prepared a GT-seq library to simultaneously genotype all 591 University and Génome Québec Innovation Centre. Raw sequencing data were genotyped using the GT-seq pipeline available on GitHub (https://github.com/GTseq /GTseq -Pipeline). The GTseq_Genotyper_ v3.pl script was used to call genotypes for each individual sample, and then, GTseq_GenoCompile_v3.pl was used to compile all the individual genotypes. This resulting output file was then converted into a PED file format (Purcell et al., 2007) for downstream genetic analysis.
Prior to analysis, we quality-filtered the entire dataset of 608 samples (591 individuals + 17 repeats). Those individuals and SNP loci with >30% missing data were removed from the dataset. Loci with a minor allele frequency <0.05 and those that deviated from Hardy-Weinberg Equilibrium (HWE) in >60% of sampling units were also removed. We calculated genotyping error for repeated samples that passed data quality filters (n = 17) as the proportion of discordant genotypes across all loci successfully genotyped for both replicates. To reduce bias due to sampling potentially highly related individuals at dens, we used ML-RELATE (Kalinowski, Wagner, & Taper, 2006) to infer relationships between all individuals using default parameters. After cross-referencing with the age class data, those individuals with an inferred parent as a result of this analysis were subsequently removed from the dataset for all population genetic analyses, except those assessing patterns of relatedness.

| Population structure and migration
To infer broad-scale patterns of population structure across the northern range of the Western Rattlesnake, we conducted Bayesian clustering analyses with STRUCTURE v 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). Ten iterations were run for K = 1-8 (the number of a priori geographic regions + 2), each with a burn-in of 500,000 and a run length of 1,000,000 MCMC steps.
The resulting output was then summarized using STRUCTURE HARVESTER (Earl & vonHoldt, 2012). To infer the optimal K value, we employed a combination of the ΔK method (Evanno, Regnaut, & Goudet, 2005) and the plotting of the log probability of the data (Pritchard et al., 2000) to assess where ln Pr(X|K) plateaued (see STRUCTURE manual). Additional iterations of STRUCTURE were also conducted for each independent cluster to test for finescale genetic substructure. For the optimal values of K, results were summarized with CLUMPP (Jakobsson & Rosenberg, 2007) and then plotted with DISTRUCT v.1.1 (Rosenberg, 2004). Based on these analyses, we defined those groups of dens that either formed distinct genetic clusters or were geographically isolated as "den complexes" (Figure 1). Using this information, we calculated Weir and Cockerham's (1984) unbiased estimates of pairwise genetic differentiation (θ) among sampling regions and den complexes as implemented in GENETIX 4.05 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004).
To further examine connectivity across the range, we tested for patterns of isolation-by-distance (Wright, 1943) among all den complexes using a Mantel test (Mantel, 1967) as implemented in the R package ade4 (Dray & Dufour, 2007) using 10,000 replicates to test for significance. To conduct this test, the centroid location of each den complex was determined using the rgeos (Bivand & Rundel, 2014) package in R (R Core Team, 2019). The straight-line geographic distance between centroids was then calculated in kilometers with GeographicDistanceMatrixGenerator_v1.2.3 (http:// biodi versi tyinf ormat ics.amnh.org/open_sourc e/gdmg/downl oad. php) and correlated with estimates of pairwise θ. We then repeated this test using den-level pairwise geographic and genetic distances within each geographic region containing >2 dens to evaluate more fine-scale patterns of connectivity.
We also assessed the amount and direction of contemporary migration among den complexes with BAYESASS 3.04 (Wilson & Rannala, 2003). This Bayesian approach uses multilocus genotypes to estimate movement among populations over the past three generations. We estimated levels of migration using 10,000,000 iterations after a burn-in of 1,000,000 MCMC steps, sampling after every 100 steps. Five runs were conducted, each starting with a different random seed, and convergence was analyzed by comparing mean posterior estimates of migration. To determine whether migration between den complexes was significant, we calculated 95% credible sets by calculating the mean migration rate ±1.96x the standard deviation, as suggested in the user manual. Those migration estimates with credible sets that did not include zero were deemed significant.

| Genetic diversity and relatedness within and among den complexes
We evaluated genetic diversity within each den complex and within each den by calculating estimates of observed heterozygosity (H o ), expected heterozygosity (H E ), nucleotide diversity (π), and the proportion of variable sites (P) as implemented in ARLEQUIN (Excoffier & Lischer, 2010). Global tests for heterozygote deficit or excess were run using Fisher's exact tests in GENEPOP 4.5 (Rousset, 2008), with 10,000 dememorization steps, 100 batches, and 10,000 iterations.
Those groups with a p-value <.05 were considered to significantly deviate from the level of heterozygosity expected if they were in Hardy-Weinberg equilibrium. We also calculated inbreeding coefficients (F IS ) for each den complex with GENETIX (Belkhir et al., 2004), using 1,000 permutations to test for significance. To determine whether measures of genetic diversity were correlated with latitude, linear regressions with nucleotide diversity and the proportion of polymorphic SNP loci (P) in a den complex were conducted in R (R Core Team, 2019) for all 11 den complexes.
Additionally, pairwise estimates of Queller and Goodnight's unbiased relatedness coefficient were calculated for individuals across each geographic region as well as within and among den complexes and dens using SPAGeDI 1.5 (Hardy & Vekemans, 2002). Reference allele frequencies were determined by pooling all individuals in a given geographic sampling region. For this analysis, we added the individuals with presumed parents back into the dataset.

| GT-seq genotyping
A total of 608 samples (591 individuals + 17 repeats) were genotyped, with an average of 341,838 raw sequencing reads per sample. After quality filtering for 30% missing data, minor allele frequency > 0.05, and HWE, 505 samples and 308 loci remained in our dataset. Mean genotyping error across all repeated samples (n = 17) was 0.36%.
After repeated samples and those individuals with an inferred parent in the dataset (n = 27) were removed, a total of 461 individuals were used for subsequent population genetic analysis with an average depth of 332x (range: 23x-1113x) and 3.82% missing data (range: 0.00%-29.87%).  (Table 1). Across the identified den complexes, estimates of θ ranged from 0.0459 to 0.526 (p < .001) ( Table 2).
Within sampling regions, there was no significant pattern of F I G U R E 2 STRUCTURE bar plots displaying inferred clustering and ancestry estimation for all individuals (n = 461), showing K = 3 for iteration 1 and K = 3 for iteration 2 for analyses run among geographic regions. Iteration 3 investigated evidence for further substructure within regions, revealing K = 3, K = 2, K = 2, and K = 2, respectively (population acronyms below plots as defined in Figure 1). Each color represents a unique genetic cluster, and each vertical bar represents the proportion of ancestry to each cluster for each individual.  isolation-by-distance detected among dens within the TNI (r = .221, p = .298) or GF (r = .911, p = .167) regions; however, a significant pattern of isolation-by-distance was detected among dens in VER (r = .873, p = .03), OKS (r = .811, p < .001), and WA (r = .57, p < .001).
Estimates of recent migration found significant gene flow from the OS den complex to WA-S and from OS to WL. All other estimates of gene flow were low (m < 0.1) and not significant (Table 3).  (Table 4). There was no significant correlation detected between latitude and proportion of polymorphism (r 2 = .196, p = .10); however, nucleotide diversity was significantly and negatively correlated with latitude (r 2 = .365, p = .03) (Figure 4).

| Genetic diversity and relatedness
At the individual den level, estimates of heterozygosity tended to be higher than those calculated at the den complex level (Table 4).
Values of observed heterozygosity ranged from 0.305 (TNI-C02) to 0.603 (WA-S01). Evidence for a significant heterozygote deficit was found in the WA-N03, WA-N04, WL-05, and TNI-C02 dens, all of which also exhibited significant levels of inbreeding (Table 4). In contrast, the WA-N06, WA-S01, VER-KL03, VER-CB01, and VER-CB02 dens all displayed an excess of heterozygotes, which was not observed in their associated den complexes overall (WA-N, WA-S, VER-KL, or VER-CB). Estimates of nucleotide diversity and the proportion of polymorphism among individual dens were similar to those determined at the den complex level (Table 4).
Additionally, average pairwise relatedness generally was higher between snakes within a den than either among dens in a den complex or among all snakes within a geographic region, though this varied across the distribution. Across regions, mean pairwise relatedness ranged from −0.05 in MW to −0.005 in OKS (Table 5). TNI-E was the most related of all 11 den complexes (r QG = .781), while MW was the least related (r QG = −.047) ( Table 5). Within dens, mean pairwise relatedness varied from −0.077 (WL-04) to 0.781 (TNI-E01) ( Table 5).

| Population structure and connectivity
We found evidence for significant hierarchical genetic structure and low connectivity among Western Rattlesnakes at the northern ex- In particular, TNI and VER each represented a unique genetic cluster and exhibited a high degree of genetic differentiation from all other sampled regions ( Figure 2, Table 1). The other geographic regions (OKS, MW, and GF) together represented a third unique genetic cluster, with evidence for underlying substructure where OKS and MW were found to be genetically distinct ( Figure 2, Table 1). In the case of GF, shared ancestry was detected with both OKS and MW, yet this region was still significantly differentiated from both based on pairwise θ (Table 1). Overall, the broad-scale patterns of genetic distinctiveness identified across the distribution generally align with identified geographic disjunction between regions due to the presence of unsuitable habitat (COSEWIC, 2015). This is especially apparent between the TNI and VER regions, where there is a high and significant level of genetic differentiation ( Table 1) that corresponds to a large disconnect in the species range ( Figure 1). This separation is believed to have existed since the Holocene Climate Optimum (Hobbs, 2013), after which temperatures began to cool roughly 6,000-5,000 years before present representing the start of a period known as the "Mesothermic Interval" in BC (Walker, 2004).
At a finer scale, we found evidence for low and generally nonsignificant estimates of genetic migration, with no significant pattern of isolation-by-distance among all 11 den complexes that we defined across the distribution (Table 3, Figure 3). While male Western Rattlesnakes will move considerable distances during seasonal migrations, constraints on migration behavior and den site fidelity may be limiting connectivity between complexes of dens across the landscape (Harvey, 2015). Adult males of this species in BC have been shown to travel a maximum of ~4 km from their den of origin prior to returning to hibernate, although migration behavior is not necessarily consistent among individual snakes or populations within the same geographical region (Harvey, 2015). Migration distances vary considerably between F I G U R E 3 Mantel's correlation between pairwise straight-line distance (km) and genetic distance (θ) for all 11 den complexes (r = −.03812, p = .4989) snakes that occupy different dens, even those neighboring one another (Gomez et al., 2015). Differences in migration patterns could potentially be a result of differences in prey availability, thermoregulatory opportunity, and habitat suitability, all of which are key factors for snake persistence (Lomas, Larsen, & Bishop, 2015). In the present study, straight-line distances between den complexes defined across the landscape ranged from 15 km (within cluster) to 412 km (among cluster) (Table 2), distances well beyond that undertaken by an individual snake in a single season. Therefore, it could take several years for snakes to encounter those from other den complexes. Moreover, environmental variability among sampled regions across the range may also be decreasing the likelihood that individuals will disperse far enough to find mates from dens outside of a given den complex. Based upon migration estimates, this seems plausible; however, as we were unable to sample from all known dens across the distribution, our sampling at the targeted den complexes may have excluded potential intermediary populations that could provide connectivity.
Estimates of relatedness within and among den complexes also indicated that the behavior of these snakes may influence levels of migration and population connectivity. Snakes sampled within a given den complex were found to be more related than snakes among complexes within a geographic region. At a finer level, snakes sampled at the same den were found, in some instances, to be more related to one another than snakes across dens within a given den

TA B L E 4
Estimates of mean observed and expected heterozygosity, proportion of polymorphism (P ), nucleotide diversity (π ), and inbreeding coefficients (F IS) for each den complex and dens within each den complex. * indicates p < .05, ** indicates p < .001 In contrast to the TNI region, our data provided weak evidence for distinct geographical or anthropogenic barriers to gene flow within the VER, OKS, and WA regions. Despite significant θ estimates between den complexes within these regions (Table 2), all values were low and fall within the level of "inbreeding connectivity" (θ < 0.20) as defined by Lowe and Allendorf (2010). This distinction suggests that gene flow may be occurring between den complexes in each of these regions at a level that can counteract the negative effects of inbreeding (Lowe & Allendorf, 2010). Within OKS, significant migration was detected between the OS and WL den complexes, which is consistent with the low estimate of pairwise differentiation. Significant patterns of isolation-by-distance were also detected among dens in both OKS and WA, suggesting that geographic distance may be the main factor responsible for limiting population connectivity throughout these regions.

C. oreganus in Canada and populations further south in the United
States (Stebbins, 2003), a significant level of gene flow was detected among the OKS and WA regions, specifically between the OS and WA-S den complexes. Given that the OS den complex is closer in proximity to the WA-N den complex, we expected to find that more migration would be occurring between these areas. While further investigation found a significant pattern of isolation-by-distance among dens across OS and WA-N (r = .518, p = .002), other factors may be contributing to the genetic differentiation detected between these den complexes. In addition to geographic distance between these sites, there are several high-elevation mountain peaks within the Pasayten Wilderness in Washington (USA) that are situated between OS and WA-N, which may be inhibiting snake movement. In order for gene flow to occur between these groups, snakes would have to successfully traverse over these mountains at elevations >2,000 m, well above the elevation at which this species has been shown to occur in studies across its range in North America (Ashton, 2001;Gomez et al., 2015;Loughran, Beck, & Weaver, 2015;Macartney & Gregory, 1988;Maida et al., 2018;Winton et al., 2018). In contrast, elevation remains relatively low (< 834m) across the landscape from OS to WA-S, potentially better facilitating connectivity between these two areas.

| Genetic diversity and relatedness
Across both den complexes and individual dens, identified patterns of genetic variation aligned with the predictions of the central-marginal hypothesis for populations that exist at the periphery of a species distribution. Here, we saw reduced levels of genetic variation within the northernmost region of the Western Rattlesnake distribution, as evidenced by a significant negative relationship detected between latitude and within den complex estimates of nucleotide diversity (Figure 4). The proportion of polymorphism detected within den complexes also tended to be qualitatively lower at the northern extent of the species range in TNI and VER, and higher in the southern populations in OKS and WA (  (Wahlund, 1928), which results from the pooling of individuals from genetically structured groups, or in this case, dens. Further analysis of dens within this group revealed low, but significant, pairwise differentiation (θ = 0.0763; p < .001) between the TNI-W01 and TNI-W02 dens, suggesting that underlying substructure may be responsible for the reduced heterozygosity of the den complex as a whole. As for the WA-S den complex, the only significant (p < .05) pairwise differentiation detected was between the WA-S02 and WA-S03 dens (θ = 0.0813). Rather than genetic substructure influencing heterozygosity at the den complex level in WA-S, relatedness among sampled individuals may be leading to a "Family Wahlund Effect" (Castric, Bernatchez, Belkhir, & Bonhomme, 2002). In this den complex, mean pairwise relatedness was higher within den than among dens, suggesting that individual dens may be comprised of separate family groups. In fact, den WA-03 contained three pairs of individuals with inferred relatedness at the level of half-siblings (data not shown). As a result, the probability of having sampled sets of closely related individuals across this entire den complex is higher than expected if individuals were randomly mating across F I G U R E 4 Correlation between latitude and nucleotide diversity (π) (r 2 = .365, p = .03), and latitude and proportion of polymorphism (P) (r 2 = .196, p = .1) across all 11 den complexes TA B L E 5 Estimates of Queller and Goodnight's (1989) pairwise relatedness (r QG ) between pairs of individuals within each region, den complex, and each den

| Implications for conservation
Taken together, our results provide valuable information for conservation and management of Western Rattlesnakes at their northern range limits, and for communal denning vipers elsewhere.
Evidence of significant population differentiation and limited connectivity within and among the main geographic regions of Western can be used to quantify adaptive genetic variation and better address the evolutionary significance criterion (Funk et al., 2012).
In addition to informing conservation unit designation, the results of this study can help guide management strategies for this and other species with similar life-history traits. For snake populations that occur in fragmented habitat, re-establishing connectivity may be challenging. Based upon the inferred patterns of hierarchical population structure and migration, the rehabilitation of fragmented landscapes at key locations may be the most effective approach to restoring population connectivity at local (i.e., the den complex level) and regional scales. The identification and subsequent preservation/ restoration of key movement corridors and critical habitats can help to encourage connectivity, whether it facilitates individual movements or the gradual expansion of existing populations. However, this approach would be challenging in landscapes where land ownership is primarily private and where strong legislation does not regulate land use to favor species at risk conservation and recovery.
Furthermore, increased public outreach and education will be necessary to minimize human-snake conflicts if these animals are to re-establish in areas from which they have been extirpated due to development.
An alternative approach for managing populations separated by unsuitable habitat or more distinct barriers could be to artificially increase gene flow through translocation. For example, in the north and central OKS region near the City of Kelowna (Figure 1), we were unable to sample snakes because population densities were too low or absent at potential denning sites. This is likely due to the intensity of human activities and development that is preventing gene flow among snakes in this area. In such cases where barriers to natural movement will be exceedingly difficult to remove, translocation could help introduce novel genetic variation across these populations and provide greater resilience to stochasticity and ongoing anthropogenic development. Translocations have been shown to be successful for both Antiguan Racers (Alsophis antiguae; Daltry et al., 2017) and Adders (Vipera berus; Madsen, Ujvari, & Olsson, 2004), resulting in bolstered population sizes and increased genetic variation (Madsen et al., 2004). However, it is important to note that in most instances, translocations are ineffective, especially for snakes and other reptiles (Brown, Bishop, & Brooks, 2009;Germano & Bishop, 2009;Nowak, Hare, & McNally, 2002;Reinert & Rupert, 1999). Many studies have found translocation to alter snake behavior and movement patterns, while also increasing population mortality (Reinert, 1991;Sullivan, Nowak, & Kwiatkowski, 2015). Additionally, the movement of individuals may also result in reduced genetic diversity for both newly established and source populations (Furlan et al., 2020), as well as outbreeding depression (Weeks et al., 2011). Due to the risks associated with translocations, it is imperative that they are biologically informed to maximize potential success. Several studies discuss the impacts of translocation on both population genetics and demography, and outline guidelines relevant to their use as a conservation tool (Furlan et al., 2020;Weeks et al., 2011). Notably, translocation efforts must account for species evolutionary trajectory as well as life-history characteristics in order to help meet conservation goals.
Overall, the continued genetic monitoring of Western Rattlesnakes in this system can allow conservation managers to track changes in genetic diversity and population connectivity over time, as well as evaluate the outcomes of management efforts.

TA B L E 5 (Continued)
The genotyping method employed in this study (GT-seq, Campbell et al., 2015) is a reliable, quick, and relatively inexpensive approach for simultaneously genotyping hundreds to thousands of individuals at hundreds of SNPs, and can be successfully applied to minimally invasive DNA samples . As the environment is continually modified due to anthropogenic development and changing climates, leading-edge genomic tools can allow wildlife and conservation managers to more effectively survey natural populations and adapt conservation strategies as necessary to encourage the persistence of species amidst their current decline.

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
Individual metadata and SNP genotypic data collected via GT-seq are deposited in DRYAD (https://doi.org/10.5061/dryad.fbg79 cns4).
All Illumina raw reads are available from the NCBI sequence read archive (BioProject ID: PRJNA631230).