Genetic signatures of small effective population sizes and demographic declines in an endangered rattlesnake, Sistrurus catenatus

Abstract Endangered species that exist in small isolated populations are at elevated risk of losing adaptive variation due to genetic drift. Analyses that estimate short‐term effective population sizes, characterize historical demographic processes, and project the trajectory of genetic variation into the future are useful for predicting how levels of genetic diversity may change. Here, we use data from two independent types of genetic markers (single nucleotide polymorphisms [SNPs] and microsatellites) to evaluate genetic diversity in 17 populations spanning the geographic range of the endangered eastern massasauga rattlesnake (Sistrurus catenatus). First, we use SNP data to confirm previous reports that these populations exhibit high levels of genetic structure (overall Fst = 0.25). Second, we show that most populations have contemporary Ne estimates <50. Heterozygosity–fitness correlations in these populations provided no evidence for a genetic cost to living in small populations, though these tests may lack power. Third, model‐based demographic analyses of individual populations indicate that all have experienced declines, with the onset of many of these declines occurring over timescales consistent with anthropogenic impacts (<200 years). Finally, forward simulations of the expected loss of variation in relatively large (Ne = 50) and small (Ne = 10) populations indicate they will lose a substantial amount of their current standing neutral variation (63% and 99%, respectively) over the next 100 years. Our results argue that drift has a significant and increasing impact on levels of genetic variation in isolated populations of this snake, and efforts to assess and mitigate associated impacts on adaptive variation should be components of the management of this endangered reptile.


| INTRODUC TI ON
Small, recently isolated populations face the risk of reduced viability over time due to genetic costs of inbreeding (Allendorf, Luikart, & Aitken, 2013). Habitat fragmentation can reduce gene flow between populations and decrease population size. This in turn increases the magnitude of genetic drift and the negative impact of inbreeding depression on survival and reproduction, and can lead to populations entering an "extinction vortex" (Gilpin & Soulé, 1986). Genetic analyses of historical and contemporary features of populations can provide information that is useful in assessing whether small isolated populations are at risk of high levels of drift and inbreeding depression, and thus should be significant components of conservation plans for endangered species (Jamieson & Allendorf, 2012).
One population feature that is important to understand in this context is effective population size (Ne) (Luikart, Ryman, Tallmon, Schwartz, & Allendorf, 2010), which measures the strength of genetic drift in a population. This is expected to reflect the rate at which variation is lost from populations and, in turn, may predict their ability to respond to future environmental change (Allendorf et al., 2013). Effective size is estimated over historical or contemporary timescales (Hare et al., 2011). Historical (long-term) estimates of Ne evaluate the parameter over evolutionary timescales and generally estimate effective size at a species-wide level (Charlesworth & Willis, 2009). In contrast, contemporary (current) measures of Ne estimate the parameter over short timescales (<5 generations) and so are most relevant to conservation and wildlife management applications (Hare et al., 2011).
The development of estimators of contemporary Ne using genetic data is an active area of research in conservation genetics (Waples, 2016;Waples, Larson, & Waples, 2016). Recently, Gilbert and Whitlock (2015) simulated data under a variety of demographic scenarios to evaluate different estimators of contemporary Ne and concluded that an estimator based on genome-wide estimates of linkage disequilibrium (LDNe, Waples & Do, 2008), as implemented in the program NeEstimator (Do, Waples, & R.S., Peel, D., Macbeth, G.M., Tillett, B.J., and Ovenden, J.R., 2014) was accurate in small isolated populations. Likewise, Wang (2016) validated the use of a sibship-based estimator implemented in the program Colony (Jones & Wang, 2010) under similar demographic scenarios. Both methods have been widely used for estimating contemporary Ne for captive and wild populations (Husemann, Zachos, Paxton, & Habel, 2016).

Analysis of observed distributions of genetic polymorphism
can also yield insights into the demographic history of populations (Lawton-Rauh, 2008;Marjoram & Tavare, 2006). To date, these analyses have primarily made qualitative inferences of changes in population size over time by comparing observed and expected patterns of diversity at small numbers of microsatellite loci (Garza & Williamson, 2001;Peery et al., 2012). However, these approaches can have low statistical power to detect even severe bottlenecks in natural populations (Busch, Waser, & DeWoody, 2007;Peery et al., 2012). Increasingly, alternative coalescent-based modeling approaches are being used (i.e., Storz & Beaumont, 2002;Excoffier et al., 2013), which have the advantage of using an explicit statistical framework to test diverse demographic scenarios. These approaches have been used to detect significant demographic events over recent timescales, including those within the range of significant human impacts to the landscape (Goossens et al., 2006;Harris et al., 2016;Hoffman, Grant, Forcada, & Phillips, 2011). When combined with other information, the inferred timing of demographic events can be used to assess whether humans have played a significant role in those events (Goossens et al., 2006) or whether they are instead more likely due to natural large-scale events, such as those associated with historical changes in climate (Tucker, Schwartz, Truex, Pilgrim, & Allendorf, 2012).
Finally, projecting future changes in genetic diversity can identify genetic risks faced by populations that have recently declined (Amos & Balmford, 2001). This issue is especially important when recent declines have led to populations not in genetic equilibrium.
Populations with a small contemporary effective size may retain the genetic profile of a larger population after the decline, and this can lead to overestimates of the amount of variation present relative to present and future population sizes. Demographic simulations can help assess the rate at which existing levels of variation will decline and provide insights into an effective timeline for management activities such as genetic rescue (Frankham, 2015;Tallmon, Luikart, & Waples, 2004). One difficulty is that while the ultimate goal of management activities is the conservation of adaptive variation, these simulations are based on neutral evolutionary processes like genetic drift and focus on patterns of change in neutral variation (Carvajal-Rodríguez, 2010). Nonetheless, they provide a time frame in which adaptive variation may be lost under the assumption that levels of adaptive and neutral variation covary (Reed & Frankham, 2003).
The eastern massasauga rattlesnake (Sistrurus catenatus) is a small short-lived (generation time ~2 years) venomous snake found throughout eastern North America. In the United States, it primarily exists in relatively small, isolated patches of habitat surrounded by heavily modified landscapes (Szymanski et al., 2016). Population declines throughout the range due to habitat fragmentation and destruction have led to the listing of this species under the US Endangered Species Act (US Fish & Wildlife Service, 2016) in the United States and as a Species at Risk (Government of Canada, 2009) in Canada. This species exhibits little phylogeographic structure across its range (Sovic, Fries, & Gibbs, 2016), and so important management units within this species may best align with existing individual populations. Previous work at the population level (Chiucchi & Gibbs, 2010;Gibbs & Chiucchi, 2012) identified high levels of genetic structure among populations that have persisted over long periods of time, large variation in long-term Ne between populations, and little evidence for recent declines in population sizes. However, these results were based on microsatellite-based variation alone and used qualitative methods of analyses for inferring demographic patterns.
In particular, measures of effective population size were historical long-term measures. These may not reflect the strength at which drift is currently operating in these populations, making them uninformative for present-day management decisions for this species.
Here, we use two independent genetic datasets to make conservation-relevant inferences about genetic and demographic characteristics of individual populations across the range of this species.
Our specific goals are to (a) evaluate the degree of genetic isolation among populations using a new genome-scale dataset; (b) estimate contemporary effective population sizes to generate a better understanding of the present and future importance of drift in these populations; (c) re-evaluate previous results (Gibbs & Chiucchi, 2012) that the genetic cost to living in small populations is small by using new estimates of individual heterozygosity from RADSeq loci in heterozygosity-fitness correlation (HFC) analyses (Chapman, Nakagawa, Coltman, Slate, & Sheldon, 2009); (d) test for evidence of past population size changes and estimate the timescale on which inferred changes occurred; and (e) use simulations to project the rate at which existing levels of variation will be lost in populations. Our study represents the most comprehensive analyses of the present-day and future genetic demography of existing populations of S. catenatus and provides information useful for developing the Recovery Plan for this species (US Fish and Wildlife Service, in preparation).

| Samples
We analyzed genetic data from S. catenatus individuals from 17 populations spanning the geographic range of this species (Figure 1; Table 1). Individual populations were defined as sets of randomly collected adult samples from within a single geographic location ≤3 km 2 in area for which sample sizes were sufficiently large to conduct the analyses described below (see Table 1 for marker-specific sample sizes for individual populations). Samples from most populations were collected across multiple years and/or consisted of multiple body size classes from a single year and so represent adults from a mixture of age classes. This is relevant to interpreting estimates of contemporary effective population size (see below). These populations were the same as those analyzed by Chiucchi and Gibbs (2010) with the difference that samples from one additional population from Spring Valley, Ohio (SPVY-see Figure 1), were analyzed using RADSeq data.

| Genetic data
We analyzed individuals (n = 385) using two multilocus genetic datasets: One based on DNA microsatellite loci previously analyzed by  Chiucchi and Gibbs (2010). For RADSeq loci, we generated data following the protocols described in Sovic et al. (2016). Briefly, high-quality genomic DNA was extracted from blood or tissue samples, and double-digest RADSeq libraries (Peterson, TA B L E 1 Estimates of genetic variation for each population based on microsatellite and polymorphic RADSeq loci. Microsatellite results are based on data from 17 microsatellite loci reported in Chiucchi and Gibbs (2010)  Note. Code is the identifier used for each population on Figure 1; N = number of individuals genotyped; Ho and He are observed and expected heterozygosity, Fis is the fixation index, and AR is allelic richness which were calculated using Arlequin (Excoffier & Lischer, 2010). N loci (100%) is the number of polymorphic RAD loci scored in individuals in that population (% scored in all individuals).
Weber, Kay, Fisher, & Hoekstra, 2012) were generated with EcoRI and SbfI restriction enzymes (New England Biolabs, Ipswich, MA) and a modified version of the RADSeq protocol described in DaCosta and Sorenson (2014). For details of the library preparation methods, see the Supplemental Information section of Sovic et al. (2016).
Pooled libraries of up to 36 individuals were sequenced in single-end 50-bp runs as partial lanes on an Illumina HiSeq 2,500 sequencer. De novo locus assembly, SNP identification, and genotyping of RADSeq loci were carried out on the raw fastq data using AftrRAD version 5.0 (Sovic, Fries, & Gibbs, 2015). Unless otherwise specified, default settings were used in the AftrRAD run. Specifically, reads containing one or more bases for which the quality (Phred) score was below 20 were removed from analyses (minQual = 20). Reads that shared 90% or more sequence similarity after accounting for indels were assigned as allelic variants at a given locus as part of the de novo assembly (minIden = 90). After assembly, genotypes were called at loci for which there was a minimum of five reads in the individual (MinReads = 5). Loci with <5 reads were scored as missing data.
Within each population, we then retained for analyses only those loci that were scored in all samples. Genotypes were assigned based on read (haplotype) identity for all analyses, with the exception of site frequency spectra used for model testing in fastsimcoal, which were constructed by selecting the first SNP from each locus (see details below).

| Population structure
Repeated analyses based on DNA-based markers such as RAPDs (Lougheed, Gibbs, Prior, & Weatherhead, 2000) and microsatellites (Chiucchi & Gibbs, 2010) have shown high levels of genetic structure and low migration rates between individual populations of Sistrurus catenatus. To assess whether RADSeq data showed similar levels of genetic structure, we used hierfstat (Goudet, 2005) to estimate pairwise and overall Fst for both the RADSeq dataset newly reported here and, for comparative purposes, for the microsatellite data reported by Chiucchi and Gibbs (2010). We calculated Pearson's r correlation and associated significance levels between pairwise Fst values generated for microsatellite and RADSeq data using the rcorr function in the R base package (R Development Core Team, 2015).

| Effective population size
We assumed that each population was essentially genetically isolated due to high degrees of genetic structure between populations and low levels of genetically effective migration inferred by Chiucchi and Gibbs (2010). As such, we analyzed each set of samples from a single location as a genetically isolated population. We calculated contemporary effective population sizes using two different estimators for both the RADSeq and microsatellite datasets. First, we estimated contemporary Ne using the LDNe method (Waples & Do, 2008), as implemented in the program NeEstimator (Do et al., 2014).
This method estimates Ne based on patterns of linkage disequilibrium between loci and was shown to perform well relative to other methods when calculating Ne under scenarios of low Ne and low migration rates (Gilbert & Whitlock, 2015). We used a "two allele" minimum for each locus within each population based on the recommendations of Waples and Do (2010) relative to the sample size of individuals (≤25) in almost all our populations. Confidence intervals for Ne values were estimated using a parametric approach implemented in the program. Second, we used Colony (Jones & Wang, 2010), which uses a maximum-likelihood-based method to calculate Ne based on inferred sibship frequencies among the samples, with associated confidence intervals obtained through bootstrap resampling.

| Heterozygosity-fitness correlations
Snout-vent length (SVL) and body mass data were available for a subset (N = 91) of the individuals genotyped for RAD variation. These (ln) these data and performed a major axis regression in the R package lmodel2 (Legende, 2018) to obtain the slope of the best-fit line between these variables.

| Model-based analyses of historical demography
We used historical demographic analyses to test between models of population history. Our goal was to assess whether evidence exists for declines in population size within each population. For declining populations, we assessed whether the declines occurred over a timescale consistent with anthropogenic impacts (arbitrarily defined as <200 years-see Discussion) or were better explained by historical factors occurring over longer timescales. To avoid problems associated with model overparameterization, we analyzed the simplest models possible that captured the key features of the process (population size change) we were interested in exploring. Specifically, we tested between two simple demographic models ( This approach has advantages over qualitative analyses of demography based on microsatellite data (e.g., Kimmel et al., 1998)  A key parameter in these models is the genome-wide mutation rate.
Because we are not aware of any direct estimates for genome-wide mutation rates for snakes, we conducted a series of preliminary analyses using paired sets of models and datasets that only varied in their use of "high" and "low" mutation rates. For the "high" mutation rate, we used 2.5 × 10 −8 per site/generation as estimated for humans (Nachman & Crowell, 2000) and used in Sovic et al. (2016). For the "low" mutation rate, we used 2.1 × 10 −9 per site/generation as used for a recent study on frogs by Thomé and Carstens (2016) and used in Gibbs et al. (2018).
We then compared the model likelihoods for each analysis and found that the use of the high mutation rate consistently generated significantly better model likelihoods (results not shown). Consequently, we used the rate of 2.5 × 10 −8 mutations per site/generation in our subsequent analyses. Finally, as described by Sovic et al. (2016), we followed the approach of Lande, Engen, and Saether (2003) to convert estimates of time from generations to years, but used a lower estimate of age at first reproduction (2 years) which is based on data collected from the populations we studied (G. Lipps, unpublished data). This resulted in an estimated generation time of 2.03 years.
We stress that although we have conducted independent historical demographic analyses on each population the degree to which each population represents independent entity over historical timescales is unclear. Thus, as a whole the results of these analyses should be taken as a qualitative survey of the types of demographic histories experienced by an unknown number of historically independent populations across the range of this species.

| Potential impact of Fst-outlier loci on demographic results
Our analyses of population structure, size, and historical demography assume that the polymorphism measured represents neutral variation, yet loci showing high levels of divergence between populations (outlier loci) may represent loci under divergent selection.
F I G U R E 2 Historical demographic models analyzed for each population using fastsimcoal. Note that for the "bottleneck model," there are no restrictions regarding the direction of the population size change, as it allows for both increases or decreases in population size. Table S3 shows the results for model comparisons for each population based on AIC.
Parameter estimates for the best-fit models are given in Table S4 We explored the sensitivity of our demographic analyses to such loci in two ways. First, we used Outflank version 0.2  with default parameter settings to identify highly divergent loci across populations. This program represents a "second-generation" Fst-outlier detection program that has much lower false-positive rates, yet comparable power as compared to earlier programs. For the Outflank dataset, we allowed up to 10% missing data at each locus. Second, because characteristics of our study system (high levels of genetic structure and small numbers of loci and populations) may make detection of outliers problematic , we also conducted a sensitivity analyses in which we chose a set of three population pairs that represent a range of sizes and interpopulation distances, and then assessed the impact of deleting a portion of high Fst loci from the dataset on demographic estimates. Specifically, we used Genepop (Rousset, 2008)

| Projecting loss of genetic variation
We assessed how existing levels of variation might change in the future for a range of contemporary effective sizes given no change in current sizes. To do this, we conducted forward genetic simulations in simuPOP with overlapping generations (Peng & Kimmel, 2005 We compared our simulation results with a widely used analytical approach for estimating loss of variation based on population genetic theory (James, 1970) that has been used for estimating changes in heterozygosity in endangered species (Amos & Balmford, 2001).
Specifically, James (1970) proposed using the equation Ht = H0 (1−1/2Ne) t to describe changes in heterozygosity over time in bottlenecked populations where H0 is initial heterozygosity and Ht is heterozygosity t generations after the instantaneous decline of a large population to size Ne. A key difference between the two methods is that compared to the simulation-based approach in SimuPOP, this analytical approach assumes no overlapping generations with all individuals only reproducing in a single generation. As with our simuPOP runs, we fixed the initial H0 at 0.3 and then projected the loss of variation using the equation above for three fixed values of Ne: 50, 30, and 10.

| Genetic markers
For the RADSeq data, a mean of 1,037,766 sequence reads were assigned to each of the 263 samples (range 216,381-4,239,460). The average numbers of total non-paralogous loci and polymorphic loci per population were 52,815 and 641, respectively (Table 1). For a given population, the median read depth ranged from 48 to 61 reads per polymorphic locus, and the percentage of polymorphic loci that were scored for all sampled individuals ranged from 25% to 60% across the 17 datasets (Table 1). These values for coverage are consistent with those reported in Sovic et al. (2016), while the number of polymorphic loci used in the analyses here is smaller than for the previous study. This is because datasets in this study represent only intra-population variation, while the datasets for the previous study included both intra-and interpopulation variation.  (2010) that these populations are genetically distinct and demographically isolated from each other.

| Estimates of effective population size
Point estimates of the effective population sizes for individual populations vary depending on the type of genetic data and method used to generate the estimates (Figure 3; Table S2). A small number of populations (RADSeq: EHSP and VNGO; microsatellites: BPNP and GLAD) yielded estimates that were 1-2 orders of magnitude greater than other populations or were estimated as infinite in size-these are coded as NA in Table S2, and omitted from Figure 3. Estimates from LDNe for RADSeq data range between 2 (JENN) to 48 (BPNP) with an overall mean (±SE) of 17.3 ± 3.5. Microsatellite-based estimates of LDNe ranged from 3 (GRL-1) to 64 (JENN) with an overall mean that is slightly higher (20.6 ± 3.6 SE). with only a few populations having larger sizes.

| Heterozygosity-fitness correlations
There was no significant relationship between our estimate of individual fitness (standardized mass index−SMI) and genome-wide individual heterozygosity estimates based on RADSeq data (r 2 = 0.001; p = 0.418). We also performed regression analyses on individual populations with a minimum of five samples, but as for the full dataset, no significant positive relationships were detected (results not shown).

| Demographic modeling
Statistical comparisons of the two demographic models in Figure 2 show strong support (AIC ~1) for the model incorporating population size change as best-fitting the data for all individual populations (Table S3) Table S4), and distributions reflect estimates from 50 bootstrapped datasets for each population. Times along the x-axis are binned into intervals of 100 generations. The bin of 0-100 generations, roughly corresponding to an anthropogenic timeframe (~200 years), is indicated with red. All estimates exceeding 2,000 generations are combined into the 2,000 generation bin (Ne = 50) to 99% (Ne = 10) ( Figure 5). This demonstrates the importance of incorporating realistic demographic parameters into models used for projecting loss of variation. Overall, these results suggest that nearly all populations are not at genetic equilibrium relative to their current effective size but have excess polymorphism which they will lose due to drift at varying rates over the next 100 years. This implies that most have undergone recent declines in census size (reflected in estimates of current effective size), but these declines have yet to influence the standing genetic variation present in a given population. These estimates are much smaller than the effective sizes based on microsatellite variation reported by Gibbs and Chiucchi (2012) for many of the same populations reported here and, in some cases, for estimates of Ne based on historical demographic modeling of RADSeq-based polymorphism by Sovic et al. (2016) (and here). This is because these studies estimate Ne over different timescales using different methodologies, and so values from each set of studies are not directly comparable (Hare et al., 2011). For example, the estimates of Ne in Chiucchi and Gibbs (2010) are based on results from

| Estimation of contemporary Ne
Migrate (Beerli, 2009) which generates estimates of long-term effective size. These reflect the amount of genetic diversity (estimated as)maintained in a population at evolutionary equilibrium due to the joint effect of drift and mutation (µ) over a period of approximately 4Ne generations ( = 4Neµ). Likewise, the estimates of Ne using fastsimcoal generated by Sovic et al. (2016) and here also represent historical estimates of effective population size over many past generations. These long-term estimates have conservation value in potentially representing a measure of population size that predates human impacts. As such, they could represent a historical target for assessing current management goals (Hare et al., 2011) or for evaluating the magnitude of recent anthropogenic impacts relative to historical population sizes (Alter, Rynes, & Palumbi, 2007 were consistent with differences in the conservation status of each species. The long-term effective population size was much smaller for the representative S. catenatus population than for the S. tergeminus population, and the S. tergeminus population was best modeled as a growing population, whereas S. catenatus showed evidence of a long-term population decline. In contrast, the short-term estimates of Ne based on single sample estimators reported here measure effective population size over much more recent timescales-approximately the period of sample collection, and hence are useful for different conservation purposes (Hare et al., 2011). Specifically, in species with overlapping generations, they generate an estimate of Ne that is approximately equal to Nb * generation time, where Nb is the effective number of breeding individuals in a given year. As such, short-term estimates are more useful for estimating short-term abundance using empirically derived estimates of the ratio of Ne/Nc, where Nc is population census size (Palstra & Ruzzante, 2008), and evaluating the degree to which genetic drift may erode the adaptive genetic potential of small populations over timescales relevant to management activities (Lynch, Conery, & Burger, 1995).
Because our analyses use mixtures of age classes of individuals from each population, they may underestimate true values of contemporary Ne (Waples, Antao, & Luikart, 2014). The degree of this bias is Illinois reported here. As described by Bradke et al. (2018), these values fall within the range of LDNe estimates found for other threatened and non-threatened snakes (see Table 3 in Bradke et al. (2018)).
A general benchmark for evaluating whether species are at shortterm genetic risk comes from the widely debated "50/500 rule" first proposed by Franklin (1980) for assessing the minimum viable size for endangered species. Briefly, the rule proposed that to avoid the deleterious effects of inbreeding depression, populations should be maintained at a minimum short-term effective size of at least 50 individuals. The validity of 50 individuals as a specific numerical target has been much discussed (Frankham, Bradshaw, & Brook, 2014;Jamieson & Allendorf, 2012) and argued to be too small to achieve its purpose (Frankham et al., 2014). Regardless, our empirical finding that most populations of S. catenatus have a contemporary Ne of <50 suggests that these populations are potentially at risk of the deleterious effects of inbreeding depression (but see Wood, Yates, & Fraser, 2016) and that conservation planning for this species needs to take this possibility into account.

| Heterozygosity-fitness correlations
Our results reinforce those of Gibbs and Chiucchi (2012) who also found little evidence of genetic costs of inbreeding based on heterozygosity-fitness correlations (HFCs) (Chapman et al., 2009). We note that the datasets used in the two studies are not independent.

Many of the same individuals (and their estimates of body condition)
were used in both analyses, with the difference being in the type and number of genetic loci used to estimate individual heterozygosity.
We echo the cautions of Gibbs and Chiucchi (2012) (Pielou, 1991;Schmidt, 1938). In contrast, the older dates correspond with large-scale environmental changes related to climate that have previously been hypothesized to impact the distribution of this (Cook, 1992)  These results are at odds with previous work that showed limited evidence for population declines in these populations using different methods of analyses and microsatellite genetic data. Specifically, Chiucchi and Gibbs (2010) used allele distribution tests implemented in Bottleneck (Luikart, Allendorf, Cornuet, & Sherwin, 1998) and

| Historical demography
only found significant results supporting recent bottlenecks in three of 16 populations. We suggest that the SFS-based tests used here are more sensitive because they employ an explicit model-testing framework, use a much greater amount of genetic data, and examine the possibility of declines over a broader timescale.
Other studies using similar approaches have found evidence for either recent or historical bottlenecks in other endangered vertebrates (Dussex, Rawlence, & Robertson, 2015;Goossens et al., 2006;Salmona et al., 2012;Tucker et al., 2012;Zhu et al., 2013), but in general, the focus is on one or a few populations of each species.
Our study differs in that it examines evidence for declines across a relatively large number of independent populations, and shows that the timing of significant declines varies across populations. This suggests that characterizing species as impacted only at historical or contemporary timescales is simplistic and that a given species may experience both types of declines. Identifying populations that have experienced significant recent declines may represent a way of prioritizing certain populations over others in terms of the impacts of anthropogenic effects on the viability of specific populations.

| Projecting future levels of genetic variability
Our simulations suggest that if there are no changes in key population parameters like size and levels of migration, most Sistrurus catenatus populations will lose 20% or more of their existing variation over the next 100 years. These conclusions are based on several assumptions.
First, they assume populations will remain no larger than their currently estimated effective size with no increase in migration. This is reasonable given that trends over time are for populations to remain isolated with the same or declining numbers Szymanski et al.., 2016).
Second, our analyses are based on the evolutionary dynamics of what we assume is neutral variation and not the adaptive variation that is used to delineate adaptive conservation units within threatened taxa (Funk, McKay, Hohenlohe, & Allendorf, 2012). This is a common problem in conservation genetic studies in that the genes that underlie adaptive variation are difficult to identify. However, Reed and Frankham (2003)  In general, these findings suggest that adaptive variation is present in existing populations of Sisturus and that high levels of drift have yet to purge this variation. However, whether the rate at which it is expected to be lost parallels the results of our simulations is unclear.
Because selection generally acts to retard the impact of drift on the loss of variation within populations, we suspect that the timeline for loss of variation established by our simulations may be conservative.
These results have conservation implications. They suggest that existing populations of Sistrurus may not yet have paid the true genetic cost of living in their current size population, but that this cost could be "paid" in the near future. As such, despite their history of living in small populations over historical timescales (Chiucchi & Gibbs, 2010), they may be moving to a new equilibrium with respect to gene dynamics, with issues to do with losses of adaptive genetic variation becoming more prominent in the near future. In other words, many of these populations may be poised to enter the extinction vortex (Gilpin & Soulé, 1986) but have yet to experience the fitness costs that will act in a positive feedback loop to drive populations to even smaller sizes. The simulations may also provide a timeline that could guide the timing of introductions associated with genetic rescue (Tallmon et al., 2004) in the event it is adopted as a viable management option . Specifically, our simulations suggest that as a guide to preserve 50% of the existing variation in the very smallest populations (Ne = 10), new genetic variation would need to be introduced within the next 40 years, whereas for the largest populations (Ne = 50), this could be delayed for as long as >100 years. Genetic introductions have been used to increase population viability in other highly inbred snake populations through the introduction and subsequent removal of radio-tagged males (Madsen, Shine, Olsson, & Wittzell, 1999;Madsen, Ujvari, & Olsson, 2004). We feel that a similar approach could be used to supplement variation in populations of Sistrurus, especially given the fact that this species represents a taxon where outbreeding depression is not expected to be a significant issue (Frankham, 2015).

ACK N OWLED G EM ENTS
We thank all those individuals who have provided samples or assisted with collections across the range of S. catenatus over the past 25 years-this work would not be possible without their help. University and the Ohio Division of Wildlife.

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

DATA ACCE SS I B I LIT Y
Raw fastq data files used in the genetic analyses are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad. jg81r80.