Intact landscape promotes gene flow and low genetic structuring in the threatened Eastern Massasauga Rattlesnake

Abstract Genetic structuring of wild populations is dependent on environmental, ecological, and life‐history factors. The specific role environmental context plays in genetic structuring is important to conservation practitioners working with rare species across areas with varying degrees of fragmentation. We investigated fine‐scale genetic patterns of the federally threatened Eastern Massasauga Rattlesnake (Sistrurus catenatus) on a relatively undisturbed island in northern Michigan, USA. This species often persists in habitat islands throughout much of its distribution due to extensive habitat loss and distance‐limited dispersal. We found that the entire island population exhibited weak genetic structuring with spatially segregated variation in effective migration and genetic diversity. The low level of genetic structuring contrasts with previous studies in the southern part of the species’ range at comparable fine scales (~7 km), in which much higher levels of structuring were documented. The island population's genetic structuring more closely resembles that of populations from Ontario, Canada, that occupy similarly intact habitats. Intrapopulation variation in effective migration and genetic diversity likely corresponds to the presence of large inland lakes acting as barriers and more human activity in the southern portion of the island. The observed genetic structuring in this intact landscape suggests that the Eastern Massasauga is capable of sufficient interpatch movements to reduce overall genetic structuring and colonize new habitats. Landscape mosaics with multiple habitat patches and localized barriers (e.g., large water bodies or roads) will promote gene flow and natural colonization for this declining species.

netic patterns of the federally threatened Eastern Massasauga Rattlesnake (Sistrurus catenatus) on a relatively undisturbed island in northern Michigan, USA. This species often persists in habitat islands throughout much of its distribution due to extensive habitat loss and distance-limited dispersal. We found that the entire island population exhibited weak genetic structuring with spatially segregated variation in effective migration and genetic diversity. The low level of genetic structuring contrasts with previous studies in the southern part of the species' range at comparable fine scales (~7 km), in which much higher levels of structuring were documented. The island population's genetic structuring more closely resembles that of populations from Ontario, Canada, that occupy similarly intact habitats. Intrapopulation variation in effective migration and genetic diversity likely corresponds to the presence of large inland lakes acting as barriers and more human activity in the southern portion of the island. The observed genetic structuring in this intact landscape suggests that the Eastern Massasauga is capable of sufficient interpatch movements to reduce overall genetic structuring and colonize new habitats. Landscape mosaics with multiple habitat patches and localized barriers (e.g., large water bodies or roads) will promote gene flow and natural colonization for this declining species.

K E Y W O R D S
dispersal, fragmentation, Island, reptile, Snake, spatial genetics, species distribution modeling ecology, and life-history characteristics affects gene flow among populations at multiple geographic scales, thus defining the functional connectivity of populations across landscapes. Even for a single species, patterns of gene flow can vary between locales, resulting, for example, in different levels of genetic structuring between populations (Coulon et al., 2010;Dudaniec et al., 2012;Jørgensen et al., 2005). We can assume that in many circumstances the historic landscapes in which a species resides have been differentially modified following European colonization, with some sites retaining more of their historic character (Caplat et al., 2016;Landguth et al., 2010;Vera-Escalona et al., 2015). Understanding the factors that influence functional connectivity among sites that retain more or less of their historic character can help us maintain historic levels of gene flow by guiding habitat management efforts for a species.
Different environmental contexts may produce varying degrees of genetic structure, such as an undisturbed site compared to a fragmented landscape (Moore et al., 2011). Although fragmentation can initially increase population densities due to crowding in remnant patches, over time populations tend to decline due to reduced habitat availability, increased edge effects, and limited migration (Fletcher et al., 2018;Haddad et al., 2017), thereby hastening differentiation among demes (Delaney et al., 2010;Yamamoto et al., 2019). The reduction of genetic variability in small, isolated populations, due to drift and lack of gene flow, can limit their adaptive potential (Frankham, 1995;Wade et al., 2017) and increase inbreeding depression, making them vulnerable to environmental and demographic stochasticity (Fagan & Holmes, 2006;Soulé et al., 1986). These effects are more prominent for less mobile species (e.g., Amos et al., 2014), thus enhancing structural connectivity may be particularly important for gene flow of less mobile organisms within fragmented landscapes; otherwise, more direct forms of intervention may be necessary (Frankham, 2015).
Throughout most of its Great Lakes distribution in North America spanning from Iowa to New York, the Eastern Massasauga (Sistrurus catenatus; Figure 1) exemplifies the ramifications of population isolation for a dispersal-limited species. Having variable, but typically small, home range sizes (25 ha (Weatherhead & Prior, 1992), 4.02 ha (Marshall et al., 2006), 1.3 ha (Moore & Gillingham, 2006), and 0.98 ha (Reinert & Kodrich, 1982)) within their wetland-associated grassland habitats has enabled Eastern Massasaugas to persist in heavily fragmented environments, particularly in the southern part of its range (Iowa, Illinois, Indiana, Ohio, and Pennsylvania), but not without genetic consequences (see below). However, similar to other ambush predators (Shine & Fitzgerald, 1996;Webb & Shine, 1998), and given their small size relative to most rattlesnake species, Eastern Massasaugas exhibit limited dispersal and movement (maximum range lengths 1-2 km; DeGregorio et al., 2011;Durbian et al., 2008) that reduces their ability to colonize new habitats (McCluskey et al., 2018).

Eastern Massasaugas are listed as Threatened in the United
States under the Endangered Species Act and in Canada under the Species at Risk Act; therefore, concerns about inbreeding, functional connectivity, and population viability have prompted several genetic investigations across the species' range DiLeo et al., 2013;Martin et al., 2021;Sovic et al., 2019). The genetic consequences of population isolation and subsequent reduction or cessation of gene flow are evidenced by low effective population sizes (<50 individuals), estimated for multiple populations throughout the range, and an expected reduction of existing population genetic variation of at least 20% in the next 100 years (Baker et al., 2018;Bradke et al., 2018;Martin et al., 2021;Sovic et al., 2019).
In regional analyses of Eastern Massasauga genetic structure from Illinois, Ohio, and Pennsylvania, Chiucchi and Gibbs (2010) found highly structured populations, even when populations were separated by relatively small distances (<7-25 km). Sovic et al. (2019) re-evaluated the Chiucchi and Gibbs (2010) data and confirmed that isolation coupled with small population size caused declines in genetic diversity for many of these populations. These results suggest a populations-as-islands model of conservation management for the Eastern Massasauga in the highly fragmented portion of its range.
Conversely, DiLeo et al. (2013) showed that genetic isolation is not ubiquitous for this species in their study of Eastern Massasauga populations along eastern Lake Huron (Bruce Peninsula and East Georgian Bay) in Ontario. They identified multiple population clusters along East Georgian Bay but cluster membership for two of these extended beyond 25 km, exceeding the scale of the structured regional populations from Chiucchi and Gibbs (2010 with few barriers (Caizergues et al., 2003;Gibbs, 1998). The Ontario populations described by DiLeo et al. (2013) are in remote areas with minimal development or barriers (though Bruce Peninsula has substantial tourism traffic and an associated road network; Reed & McKenzie, 2010); therefore, Eastern Massasauga movement should be less restricted on Bruce Peninsula compared to the fragmented landscapes sampled by Chiucchi and Gibbs (2010) for their regional analyses.
The question remains as to whether population isolation and strong genetic structuring are the norm for Eastern Massasaugas, or whether relatively intact landscapes do broadly promote connectivity. If the latter, targeted habitat management in fragmented landscapes would have the potential to restore functional connectivity and enhance population viability. Here, we determine whether gene flow and low genetic structuring are present elsewhere in the Eastern Massasauga's range. Our study site, Bois Blanc Island (BBI), is located in Lake Huron, between Michigan's Upper and Lower Peninsulas. BBI resembles the part of Ontario with the highest abundance of Eastern Massasaugas, with numerous occupied discrete habitat patches, minimal development aside from roads, and is situated at the species' northern distribution limit (Szymanski et al., 2016). Our spatial sampling is comparable to the fine-scale ing the capacity of Eastern Massasauga to maintain gene flow over larger areas but did not evaluate spatial genetic patterns pertaining to structuring at a scale representative of interpatch movement for this species. Therefore, our study aimed to fill this gap in knowledge and represents the first effort to examine fine-scale spatial genetic structure for Eastern Massasaugas inhabiting a landscape with high patch occupancy and abundance.
Specifically, we aimed to (a) quantify available habitat and landscape connectivity across this comparatively unmodified landscape and (b) assess spatial genetic relationships using approaches not employed in previous Eastern Massasauga studies, yet at a spatial scale where strong genetic structuring has been found in Illinois, Ohio, and Pennsylvania .
We incorporated species distribution modeling (SDM) to reveal the extent and connectedness of habitat on the island so our genetic results could be evaluated within the context of habitat availability. An unmodified landscape by itself may not be sufficient for promoting successful gene flow in a dispersal-limited species without a robust habitat network. Given the similar environmental contexts (i.e., low human population density, contiguous tracts of natural vegetation), we predicted Eastern Massasauga populations on BBI would exhibit limited genetic structuring, similar to the Lake Huron region, Ontario populations. Developing Eastern Massasauga management strategies that create or maintain connectivity to promote migration and natural colonization is a priority. Our study assessed genetic processes within a landscape with high Eastern Massasauga abundance, low human development, and lacking active management, thus offering insights into Eastern Massasauga dispersal patterns in landscapes with minimal human intervention.

| Study area
Bois Blanc Island (BBI) is one of the larger islands (88 km 2 ) in the Great Lakes, situated in Lake Huron between the Upper and Lower Peninsulas of Michigan, USA. (Figure 2a). Predominantly forested with wetlands and beaver meadows interspersed throughout, the island is largely undeveloped with a small resident human population (95 people; 0.75 people km −2 (U.S. Census Bureau, 2010). The island is estimated to have been deglaciated 11,200-11,000 years ago (Larsen, 1987), after which BBI was isolated from the Upper and Lower Peninsula until 10,300 years ago when lake levels dropped, connecting BBI to the Lower Peninsula (Larsen, 1987).
This period is likely when Eastern Massasaugas colonized BBI.
The water level of the lakes once again rose and BBI became isolated around 4,500-4,000 years ago (Larsen, 1987 1978, 1979, 1990, and 2010 and by identifying open-canopy areas using aerial imagery (ArcMap Basemap; Google Earth). We distributed survey sites to ensure that the heterogeneity of the landscape and genetic variation of individuals was adequately sampled (Figure 2b; Balkenhol et al., 2015). At the location of each snake capture, we recorded soil temperature, shaded air temperature, and cloud cover. We measured each individual's weight and length and probed the cloaca and palpated snakes to determine sex and reproductive status. Individual snakes were marked with a subdermal passive integrated transponder (PIT) tag.

| Sampling methods
All snakes were scanned prior to processing for the presence of a PIT tag, to avoid resampling the same individual. Up to 200 µl of blood was collected from the caudal artery of each individual and stored in 95% ethanol. After completing data collection, snakes were returned to their capture site and released. All equipment that contacted a snake was either sterilized with a 10% bleach solution or single-use equipment was changed out between individuals, to prevent the

| Laboratory analyses
We extracted DNA from blood samples of 102 unique Eastern Massasaugas using QIAGEN DNeasy Tissue kits following the manufacturer's protocol. We genotyped individuals at 16 microsatellite loci developed by Anderson et al. (2010). Each 10 µl PCR reaction Scu209 (62°C). To detect any potential contamination of our PCR runs, we used a negative control for each reaction set. Following PCR amplification, fragments were analyzed on an ABI3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA), at Annis Water Resources Institute, Grand Valley State University. We manually scored fragments using PeakScanner (vers. 2.0). Of the 102 samples, we reamplified and genotyped ~10% (11 individuals) of our total sample to verify our genotyping results and to calculate a PCR and allele scoring error rate.

| Genetic analyses
We tested all samples for departures from Hardy-Weinberg equilibrium and linkage disequilibrium using exact tests in the program GENEPOP (vers. 4.2) (Rousset, 2008). Using the program GenAlEx (vers. 6.5) (Peakall & Smouse, 2005), we calculated the number of alleles, effective number of alleles, observed heterozygosity, and expected heterozygosity. When investigating spatial genetic patterns, it is important to interpret results in light of the isolation-by-distance (IBD) pattern that is present (Perez et al., 2018). To assess patterns of isolation by distance (IBD), we used Mantel tests within the R (vers. 3.6.1) package ADEGENET (Jombart, 2008).
We tested for spatial genetic patterns using spatial principal component analysis (sPCA) implemented in ADEGENET (Jombart, 2008). e chose the Gabriel graph option for the spatial network as it represents a compromise between Delaunay triangulation that includes more connections than is likely for Eastern Massasauga dispersal and overly restrictive nearest neighbor options (Teich et al., 2014).
Unlike Bayesian clustering programs such as STRUCTURE (Falush et al., 2003;Pritchard et al., 2000), sPCA does not assume Hardy-Weinberg equilibrium and linkage equilibrium. sPCA uses Moran's I to identify patterns of spatial autocorrelation and is effective at detecting cryptic genetic structure (Jombart, 2008;Schwartz & McKelvey, 2008;Vergara et al., 2015). The sPCA components can reveal both global (positive eigenvalues) and local (negative eigenvalues) structures. When global structure is significant, it can indicate either clusters or clines in the dataset compared to between-individual genetic differences reflected by the local scores.
We assessed both patterns with the "spca_randtest" function recommended to increase statistical power (Montano & Jombart, 2017) using 9,999 permutations.
We used two methods to identify the most likely number of genetic clusters among samples on BBI, first with the Bayesian clustering program STRUCTURE (vers. 3.4.1) (Pritchard et al., 2000).
Our STRUCTURE analysis included 20 independent runs for K = 1 through K = 15 (approximate number of sampling areas) with an initial burn-in of 50,000 followed by 500,000 iterations under the admixture model without assigned population information. We determined the most likely number of clusters using web-based StructureSelector (Li & Liu, 2018) Figure A1). We also used a multivariate method, discriminant analysis of principal components (DAPC) implemented in the R package ADEGENET (Jombart, 2008). Like sPCA, DAPC does not rely on traditional population genetic assumptions. DAPC is effective at discerning genetic clustering in more complex population genetics models (Jombart et al., 2010).
We also evaluated spatial genetic patterns on BBI with EEMS (estimated effective migration surfaces). EEMS is a visualization tool, developed for use in systems where IBD is present but deviations due to barriers or corridors have also contributed to rates of historical gene flow. Effective migration rate is based on a steppingstone model under idealized settings that would generate the same genetic dissimilarities between demes as seen in the data (Petkova et al., 2016). The effective diversity rate assesses genetic dissimilarities that exist between individuals within a single deme (Petkova et al., 2016). The EEMS program maps effective migration surfaces where genetic similarity decays more quickly than expected under IBD, indicating reduced gene flow (Petkova et al., 2016;Silliman, 2019). Similarly, effective diversity surfaces indicate regions where genetic diversity is higher or lower than average. We calculated descriptive statistics for the high and low genetic diversity regions from this analysis using the R packages STRATAG (Archer et al., 2017) and PopGenReport (Adamack & Gruber, 2014).
Our EEMS analysis was focused on the eastern two-thirds of the island. We first adjusted most default parameters so acceptance rates fell in the recommended 10%-40% range. We left the negBi-Prob parameter at the default value to maintain more tiles for each run. Then, we ran three independent chains for 40, 60, 80, and 100 demes for 2,000,000 MCMC iterations with 1,000,000 iterations of burn-in, thinning every 10,000 iterations followed by visualization with the included rEEMSplot R plotting package (https://github. com/dipet kov/eems).
We tested for spatial autocorrelation among the 102 samples to identify patterns of fine-scale genetic structuring. Using the program GenAlEx, we used a pairwise matrix of genetic and geographic distance to calculate a correlation coefficient (r) for each distance bin. Distance bins were set at 0.5 km and extended from 0 to 5 km.
Significant spatial autocorrelation occurred when r was greater or less than the 95% confidence intervals (Peakall et al., 2003).

| Habitat analyses
We developed a BBI-specific SDM for habitat availability and connectivity estimates to complement our genetic results. Our SDM was based on an ensemble of small models (ESM) approach implemented in the R package ecospat (Breiner et al., 2015;Cola et al., 2017) that was developed for situations where few occurrences exist. ESMs consist of independent models built with each pair of variables before creating an ensemble of all models with user-specific thresholds or weighting options. We created our occurrence dataset (n = 28) by spatially thinning all BBI Eastern Massasauga records obtained from the Michigan Natural Heritage Database (1978, 1979, 1990, and 2010) maintained by MNFI and our surveys using a minimum distance of 200 m to avoid overemphasizing environmental data from any individual site. Background point selection was limited to a 10 km buffered region around all occurrence data (clipped to the extent of BBI). We screened a range of candidate environmental predictor variables representing land cover, elevation, and hydrology features potentially important for Eastern Massasauga presence (Supplement 1). First, we used the "corSelect" function from the fuzzySim R package (Barbosa, 2015) to remove highly correlated variables (>0.75), retaining the higher scoring variable from each correlated pair. We further reduced variables using the jackknife of variable importance and training gain from Maxent. These variables were retained for the ESMs: 1 km topographic position index, distance to high wetland potential (merged classes (5-8) from Coastal Change Analysis Program (C-CAP) wetland potential 2017 (NOAA)), percentage of emergent wetlands (Michigan C-CAP land cover 2016 (NOAA)) within 300 m, canopy cover (Coulston et al., 2013;Yang et al., 2018) represented by four equal intervals, standard deviation of canopy cover (Coulston et al., 2013;Yang et al., 2018)  We generated ESMs using Maxent (Phillips et al., 2006;Phillips & Dudík, 2008) and artificial neural network (ANN) model types separately, then an ensemble combining both. Both Maxent and ANN have been shown to work well in an ESM framework (Breiner et al., 2018). We applied internal ecospat tuning for both Maxent and ANN. All models (Maxent, ANN, ensemble) were run 10 times with fivefold cross-validation and weighted using Somers' D (Breiner et al., 2015). We evaluated the three ESM model types with the Boyce Index (Hirzel et al., 2006) that was calculated using the "ecospat.boyce" function for the continuous output of each model within the 10 km background extent. We also compared the three ESM models to three broader SDMs created with the northern and statewide Michigan Eastern Massasauga occurrence data (including some BBI occurrences) to determine the best fit for the BBI data.
These comparisons were also based on Boyce Index values calculated within the BBI model background extent.
We applied the threshold that maximized the true skill statistic

| Genetic analyses
Across our 102 samples, Scu209 was removed from our analyses because it was monomorphic. Using the 15 remaining loci, we calculated an allele scoring error rate of 2.98%. After simple Bonferroni correction, all loci were in Hardy-Weinberg equilibrium, while 8 of 105 pairs of loci showed significant linkage disequilibrium. The number of alleles per locus ranged from 2 (Scu206) to 11 (Scu211; Scu216) ( Table 1). The average number of alleles across loci was 5.87, and the number of effective alleles was 3.0. Observed heterozygosity across loci ranged from 0.17 (Scu208) to 0.79 (Scu215) ( Table 1).
Overall, the average observed heterozygosity for BBI was 0.58 while the expected heterozygosity was 0.60 (Table 1).

We detected significant IBD for Eastern Massasaugas on BBI
(Mantel r = 0.08). The IBD results are more characteristic of a clinal pattern than defined clusters (Figure 3). Tests for both local and global structure, using sPCA, were nonsignificant ( Figure A2).
Genetic partitioning is not evident among the BBI sampling areas.
The Evanno method and DAPC selected two and three clusters, respectively. Estimates of MedMeaK, MedMedK, MaxMeaK, and MaxMedK (Puechmaille, 2016) were similar to these results, identifying K = 2 (MedMeaK; MedMedK) and 3 (MaxMeaK; MaxMedK) as the best fit for the data, with cluster membership displaying no clear spatial pattern, consistent with weak structuring and isolation by distance (Figures A3-A6). The cluster from the southeast corner of the island exhibited the most consistent assignment of individuals to one cluster.
EEMS showed that effective migration and diversity exhibited similar geographic patterns of high and low areas relative to the mean within our study area on BBI (Figure 2c,d). Effective migration exceeds what is expected under strict IBD in the northwest, and elevated genetic diversity is found in two separate sections of northern BBI. Reduced migration and diversity are concentrated in the southeastern part of BBI, between the two inland lakes and the Lake Huron shoreline (Figure 2c,d). Genetic diversity estimates from the high-and low-diversity zones (Figure 2) showed higher observed and expected heterozygosity in the high diversity areas (Tables A1-A3).
There was slightly higher allelic diversity overall in the low-diversity region though among more individuals (n = 46) than the western (n = 19) and eastern (n = 11) high diversity zones.
We identified significant positive spatial autocorrelation between individuals at approximately 1 km on BBI, suggesting Eastern Massasaugas show restricted dispersal within 1 km (Figure 4).

| Habitat analysis
Maxent models performed slightly better than ANN and the ANN/ Maxent ensemble across the 10 replicate runs (Table A4a-

| D ISCUSS I ON
Our results confirm that gene flow is possible across Eastern Massasauga populations, even at relatively broad scales (~7 km). Our well-connected landscape with abundant suitable habitat facilitated sufficient gene flow to result in very limited genetic structure, which is not typical for this species elsewhere (but see DiLeo et al., 2013).
Our study took place on an island, but our genetic diversity measures are on par with previous estimates from highly structured populations throughout the range. We show the importance of considering landscape context when inferring spatial genetic patterns, even for a species whose low vagility and movement patterns predict strong genetic structure.

| Genetic structuring
Low observed genetic structuring across our island landscape is found high levels of genetic structuring across the range. In their fine-scale regional analysis, Chiucchi and Gibbs (2010) found strong, significant genetic structuring between pairs of sites in Illinois, Ohio, and Pennsylvania at distances equivalent to, or shorter than, the extent of our BBI study area. This discrepancy of results between our analysis and Chiucchi and Gibbs' (2010) sites can likely be attributed to differences in habitat quality, quantity, and connectivity.
The landscapes surrounding the Chiucchi and Gibbs (2010) regional sites are heavily fragmented by agriculture and roads, thereby limiting overall habitat availability and restricting migration between occupied patches compared to BBI with abundant occupied habitat patches that are well aggregated and in close enough proximity (mean distance < 300 m) to facilitate interpatch movement.
The low levels of structuring found in northern populations Yu & Lu, 2018). The Great Lakes region has lost over 50% of its wetland habitat since European colonization but those losses have not been evenly distributed: Illinois (90% lost), Indiana (87%), Ohio (90%), Michigan (50%), and southern Ontario (90%) (Dahl, 1990;Snell, 1987;Suloway & Hubbell, 1994). We detected signatures of dispersal up to 1 km, which is similar to the average maximum range length measured in Eastern Massasaugas in the Bruce Peninsula, Ontario (Weatherhead & Prior, 1992), and within the range from several other studies (Szymanski, 1998). Our results indicate genetic spatial autocorrelation occurs at fine scales; therefore, the lack of strong observed genetic structuring suggests there is gene flow across the sampled area on BBI facilitated by suitable habitat acting as stepping stones (Crandall et al., 2012).

Historically, Eastern Massasauga movement would have been aided
by the decades of logging that occurred around the turn of the 20th century (Sanborn et al., 2012;Whitney, 1987) resulting in new corridors and open-canopy habitat patches. While large-scale logging has ceased, the current BBI landscape is more permeable to movement than more heavily impacted habitats found elsewhere in the species' range, which is facilitated by aggregated habitat patches (Figure 2b).
Contiguous swaths of forest are not ideal habitat, but gaps created by storms and tree falls provide thermoregulatory opportunities for snakes moving between habitat patches (Robillard & Johnson, 2015; author obs.).
While limited genetic structuring indicates an abundance of snakes in well-connected habitats, we did find some evidence of barriers to movement. Genetic structure results must be carefully scrutinized in the presence of IBD (Perez et al., 2018), and the STRUCTURE genetic clusters seem to reflect an IBD (clinal) pattern rather than discrete spatial genetic clusters ( Figures A5-A6). Eastern

| Genetic diversity
The Massasauga population on Bois Blanc Island shows a level of genetic diversity that is on par with highly structured populations throughout the rest of the range. We might expect our island population to show reduced genetic diversity, relative to mainland sites, because of the combined effects of a potential founder effect (following colonization ~10,000 years ago), no migrants from the mainland (since the island was isolated ~4,500 years ago, Larsen, 1987), and genetic drift. However, our average observed heterozygosity per locus) that are almost two times higher than ours (5.9 alleles per locus). The elevated genetic diversity seen in these Ontario populations may reflect their broad distribution and higher levels of connectivity compared to the southern part of the range.

| Management implications
The findings suggest that Eastern Massasauga gene flow is possible when suitable habitat is abundant, and habitat patches are well connected. The short nearest neighbor distance (mean = 245 m) revealed by our habitat analysis illustrates the importance of minimizing spatial separation of patches, particularly in forested landscapes.  (Martin et al., 2021;Sovic et al., 2019). As genetic drift proceeds in isolated populations, it is estimated that in the next 100 years the genetic variation of Eastern Massasauga populations will decline by 20% (Sovic et al., 2019). This trend of declining genetic variation is further supported by the loss of 67% of rare alleles over a ten-year time period from the single remaining Eastern Massasauga population in Illinois (Baker et al., 2018).

| CON CLUS IONS
Spatial genetic structure studies on snakes are still rare despite a clear need for such information, given the global decline of reptiles (Gibbons et al., 2000;Todd et al., 2010). We used multiple analytic methods to assess fine-scale genetic structuring for a species largely considered to be restricted to isolated populations due to habitat loss and limited mobility. We showed that in a landscape context with few anthropogenic disturbances, this species is capable of sufficient rates of movement and gene flow to prevent strong genetic structuring from arising at the same spatial scale at which highly structured populations occur in the southern part of the range. Restoring landscapes to resemble BBI may be difficult across much of the Eastern Massasauga's range, as this remote location has minimal human development pressures. However, our study is a demonstration that the migration ability of Eastern Massasaugas, coupled with habitat restoration that improves connectivity, could promote natural dispersal and colonization. This might be especially important for conservation in those parts of the species' range with more fragmented landscapes than BBI. There have already been genetic consequences for several of the habitat island populations in the southern part of the Eastern Massasauga's range (Martin et al., 2021;Sovic et al., 2019) and expanding the occupied area might require population augmentation measures, with considerations for local adaptation, to revitalize the gene pool in those areas to provide insurance from future stochastic events that might tend to diminish genetic variation.

ACK N OWLED G M ENTS
We thank Brendt Sheridan, Arin Thacker, Alyssa Swinehart, Jenny Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. Fish & Wildlife Service have both expressed concern over the release of any Eastern Massasauga location data.