Using Ecological Niche Models for Population and Range Estimates of a Threatened Snake Species ( Crotalus oreganus ) in Canada

: Modelling the distribution and abundance of species at risk is extremely important for their conservation and management. We used ecological niche models (ENMs) to predict the occurrence of western rattlesnakes ( Crotalus oreganus ) in British Columbia (BC), Canada. We applied this to existing population estimates to support a threshold of occurrence for management and conservation. We also identified predictors influencing rattlesnake distribution and abundance in this region. Using a Geographic Information Systems platform, we incorporated ENMs, capture–mark– recapture (CMR) and radio-telemetry results, province-wide observations, Landsat imagery and provincial databases for agricultural land use to produce quantitative, spatially explicit, population estimates across BC. Using available western rattlesnake habitat estimated at 183.9 km 2 and averaging estimates calculated from densities in three study populations, we generated a mean adult population size of 9722 (±SD 3009; 0.8 relative index of occurrence [RIO] threshold). Only a small area (21.6 km 2 ) of suitable land cover was located within protected areas, potentially protecting an estimated 1144 (±354) adults. Most suitable land cover was within 500 m of roads (170.6 km 2 ), repre-senting potential habitat being used by an estimated 9017 (±2791) adults. At the threshold RIO value chosen (0.8), only a very small area of farmland provided suitable land cover. Our results highlight the possibility of high mortality rates for western rattlesnakes near roads and the fact that protected areas do not provide sufficient coverage to conserve the population. Given that this species has relatively low mobility and high site fidelity to home ranges, our population estimate for BC provides a useful reference for the northern part of the species’ range. It also fulfills a need to estimate population size within political jurisdictions where conservation management decisions are made, as well as presenting a method that can be applied to other parts of the range, including the southern United States. Our study provides an important benchmark for future monitoring of western rattlesnakes in BC using a repeatable and transparent approach. Similar applications can be extrapolated and applied for other threatened species to identify and quantify population distributions and threats, further supporting conservation prioritization tools to be used to maximize the effectiveness of conservation strategies under financial constraints.


Introduction
Unprecedented declines and extinctions are negatively affecting global biodiversity [1][2][3][4]. These declines are likely to be exacerbated by the recent coronavirus pandemic, which is leading to increased destruction of natural land cover and illegal wildlife trade [5]. The root causes of biodiversity declines are multifaceted, with the two most apparent drivers being resource extraction and agriculture, followed by urban and industrial development, invasive species and the spread of disease [6,7]. Cumulative effects are important to consider as well; for example, climate change is predicted to alter species' distributions as well as the magnitude and severity of other threats [8,9]. Thus, it is increasingly important in systematic conservation planning [10,11] to map species and other components of biodiversity so that management options can be spatially prioritized and threats quantified and mitigated [12][13][14].
Species that are least resilient to mounting environmental perturbations are most at risk [8]. In the case of reptiles, six major threats are "habitat loss and degradation, introduced invasive species, environmental pollution, disease, unsustainable use, and global climate change" [15,16]. While no species of reptiles are endemic to Canada, all Canadian snake species reach their northern range limits within the borders of the country. These peripheral populations thus experience naturally-occurring constraints, such as non-annual reproduction, low population densities, and short activity seasons [17,18]. At-risk populations with ranges that straddle the Canada-United States border are considered especially problematic-likely because synchronizing efforts to maintain population health and manage species on each side of the border is difficult-and maintaining connectivity range-wide is key to resilience [19][20][21]. Protecting these peripheral populations is critical for three reasons: first, many of these species have relatively small ranges [21]; second, many wildlife populations in North America have a tendency to collapse towards the northern limits of their distribution [22]; and third, peripheral populations may display genetic traits not apparent elsewhere in the range of the species [23].
Conservation planning and assessment for snakes is relatively difficult [24]. Members of the suborder are often rare, secretive, cryptic and active under relatively limited environmental conditions (annually), while occurring in a wide diversity of vegetation and land cover types [25]. More than 82% of snake species worldwide lack sufficient information to even assess their conservation status [26], and long-term studies are rare (e.g., [27,28]). Such information is needed urgently because reptiles are among the most rapidly declining taxa globally [29][30][31][32][33] and yet they play key roles in ecological food webs [27]. Snake distribution is strongly associated with spatial thermal homogeneity, and everevolving anthropogenically-altered lands may facilitate drastic alterations in the heterogeneity of thermal properties on the landscape [34][35][36]. Furthermore, many snake species have a unique combination of life history characteristics such as delayed sexual maturity, high adult survival and longevity with low fecundity [37]. In addition, many species have strong spatial fidelity to, and communal use of, traditional denning sites in which to overwinter (hibernacula), as well as summer foraging grounds and seasonal migration corridors between landscape features used multiple times over an annual cycle (e.g., western rattlesnake [38][39][40]; prairie rattlesnake-Crotalus viridis [41]; eastern diamondback rattlesnake-Crotalus adamanteus [37]; garter snakes-Thamnophis atratus [27]. This renders snakes highly vulnerable to landscape perturbations or direct mortality [42]. It magnifies the need for approaches that effectively assess population size and status, and the importance of identifying critical habitat necessary for long-term persistence. Population declines can be difficult to detect and quantify; hence, long-term studies of natural populations and communities are generally regarded as indispensable for understanding normal population trends and fluctuations [43]. The challenge in demonstrating whether observed trends in estimated population sizes constitute normal fluctuations or "unnatural" declines is that most field studies of herpetofauna lack the duration or consistency to make such determinations convincingly [44]. In Canada, the western rattlesnake (Crotalus oreganus) occurs in British Columbia (BC), as the province represents the northernmost limit of the species' range which extends south through Washington, central Idaho, northern California and Oregon (Figure 1a [45]). Less than 5% of the species' North American range occurs in Canada, and within the country all five geographic region subpopulations are disjunct from each other and two (Thompson-Nicola, Vernon; see Figure 1) are now disjunct from populations in the United States [46]. Conservation of this species is further complicated by the fact that rattlesnakes in Canada are persecuted (probably more so than non-venomous snakes), and in addition, the species is cryptic and widely distributed in a variety of vegetation types and elevations [46][47][48].Although globally listed as 'Least Concern' by IUCN [45] in Canada, the western rattlesnake is listed on Schedule 1 of the federal Species at Risk Act (SARA) as a federally-threatened species at risk [46]. In BC, western rattlesnakes generally occur in grassland ecosystems which are largely restricted to valley bottoms [35,46,49] and are comprised largely of the Bunchgrass, Ponderosa Pine, and Interior Douglas-fir Biogeoclimatic Zones (BGC [50]. Western rattlesnakes have also been found in the Engelmann Spruce-Douglas-Subalpine Fir zone [51]. Only a small portion (estimated at 15%) currently lies within protected areas [52]. Despite numerous threats, relevant and accurate population estimates and contemporary demographic parameters, as well as detailed assessments of remaining western rattlesnake habitat at local, regional or range-wide extents are limited [28]. Without accurate population estimates and information on habitat availability for this imperiled species, implementing and monitoring conservation efforts across its range continue to be compromised.  Table S1 for roads). (b). Map showing current range of the western rattlesnake (Crotalus oreganus) in British Columbia, Canada, with protected areas and agricultural land overlaid. Stars show location of intensive capture-mark-recapture study areas and data sources for our models; (see main text and Table S1 for sources).
Here, we use three machine learning ecological niche models based on occurrence data and Landsat imagery to estimate the range-wide available habitat and population size of the western rattlesnake in BC. Ecological niche models (or ENMs) provide information for establishing species boundary mapping [53] and are extremely useful for monitoring the status of species at risk, as well as for systematic conservation planning efforts [54]. Although the terms species distribution modelling (SDM) and ENM are often used to describe the same concept, Feng et al. [55] suggested that the difference was that SDMs tend to focus more on the geographic distribution of species, while ENMs are more strongly geared to analyzing components of the fundamental ecological niche. These models have been used successfully to map species distributions for other snake species, including western rattlesnakes in California and Oregon, U.S.A. [56], midget faded rattlesnakes (Crotalus oreganus concolor) in Wyoming [57], eastern hog-nosed snakes (Heterodon platirhinos) in Ontario, Canada [58], invasive frogs in Hawaii [59] and amphibians in Europe [60]. Specifically, Spear et al. [57] used MaxEnt to determine factors relating to abundance of rattlesnakes, including distance to urban areas. Moreover, the importance of validating the predictions produced by ecological niche models with both empirical distributions and expert opinion was recently highlighted by Sarquis et al. [61] in their study of Bothrops alternatus, a venomous snake species in Argentina.
One challenge with using ENMs is deriving estimates of abundance and thresholds of occurrence [62]. Another is that ecologists often must choose between conducting small-scale local intensive studies over long periods to obtain detailed demographic information, versus more extensive studies using locational data over larger areas of the species' range (often lacking true absence records and with few repeated observations). Here, we address these issues by combining observational data with climate and land cover GIS layers to quantify the area of habitat that could be used by western rattlesnakes throughout BC, with density estimates from capture-mark-recapture (CMR) studies on western rattlesnake populations [28,63]. Because we know that snakes are vulnerable to the impact of roads [64] and bio-accumulation of pesticides within the study area [65] we also examined the proportion of land cover suitable for western rattlesnakes that occurs within and outside the following designations: protected areas, a road buffer and agricultural land. Our specific goals were to determine: (1) the first approximate estimate of available habitat and population size for the western rattlesnake in BC based on observational data, and use supporting estimates from den counts to determine a threshold of occurrence; (2) the effects of different land cover and topographical variables on presence of western rattlesnakes; and different scenarios-including or excluding protected areas, roads and a 500 m buffer, and agricultural land-on models of the BC population estimate. Our approach enabled us to utilize both the concepts of prior assessments from den counts, intensive studies at specific locales, and observational data from throughout the range in BC to provide the best possible large-scale population estimate for this species in Canada to date.

Western Rattlesnake Distribution in Study Area
Five putatively disjunct subpopulations of western rattlesnakes occur in BC [23], within valley systems of the southern part of the province shown as different polygons based on observational data from the BC Conservation Data Centre in Figure 1a ( Figure  1b shows protected area and agriculture overlays-see later). The three largest portions of the range fall within the watersheds of the Thompson-Nicola Rivers (in the northwest), and Vernon, together with the Okanagan-Similkameen Rivers (the main southern part of the range), in addition to two smaller populations that occur in the Midway/Rock Creek and Grand Forks area (in the southeast of the province; Figure 1a).

Species Observations
For the spatial analysis and population estimate, we used observational data gathered throughout the geographical range archived by the BC Conservation Data Centre [66] and based on the Wildlife Species Inventory [67] databases. The government of BC houses reports and spatial data resulting from research projects. These can include location points resulting from surveys (including radiotelemetry and capture-mark-recapture studies), incidental observations (e.g., citizen science-Conservation Data Centre observations confirmed by experts and Wildlife Species Inventory observations by experts/biologists, external published reports, and some direct submissions) and museum collections. As throughout their range [68][69][70][71][72], in BC, the western rattlesnake relies on hibernacula to overwinter. Following egress in the spring, the snakes migrate to, and move throughout, their summer habitat to forage and mate, returning to their hibernacula in autumn [40,49,73]. Snake locations away from dens were linked to hibernacula within a straight-line distance of 3.6-4 km radius based on previous studies [35,48], so the data we used represented observations at den sites (one den would count as one location), as well as encounters away from dens during the active season summer months. Because of restrictions on use of data from Indigenous and private land, we could only use den data from public and government federal lands. Therefore, for this study, we only used data from approximately one-quarter of the known den sites in BC. These dens were spatially dispersed throughout the rattlesnake's range in BC. Thus, both dens and observational data away from dens were included in our models.
Of the total data points obtained from the Conservation Data Centre, 1702 were survey observations (targeted surveys which could include observations where surveys were performed but no snakes found, up to 2019), 1002 were incidental (ad hoc) observations (up to 2019), and 145 were data points obtained from individual snakes monitored during radio telemetry studies [74,75]. There was a total of 85 Conservation Data Centre element occurrences that could potentially represent one or more dens. Although we used all available data (1900s to 2019), most observations (60%) were post-2002. Observations varied in accuracy. For example, point observations collected using Global Positioning Systems, GPS (post-1997) were of greater accuracy (3 m, often with a 10-12 m buffer), whereas element occurrence polygons from historical observations had lower accuracy (250 m, or in a few rare cases, 1 km). We used the centroids of these polygons as our observation points.
We also tried to ensure that all observations were as relevant as possible to contemporary land cover (rattlesnake occurrences may have declined [45] or rattlesnakes may have been extirpated from many historical locations due to continued anthropogenic development) using the following steps. First, we removed duplicate and/or overlapping observation locations using spatial analyses in R [76]. Snake observations clustered in the same area could represent records of the same individual snake (i.e., multiple counting) or individuals aggregated at den sites. To alleviate this potential bias, we used the package 'spThin' [77] to select records that were a minimum distance apart. Based on average daily movements of rattlesnakes [49] and studies of other reptile species, we selected a minimum distance of 500 m between sightings to minimize the probability of counting the same individual twice and spatial autocorrelation issues. Some researchers working on rattlesnakes or other snake species have used a larger distance (1 km) between sightings [56,78]. However, using distances > 500 m in our study would have substantially reduced sample sizes. The R package spThin deployed a randomization nearest-neighbor distance (NND) to select spatially independent observations, meaning that the filtered data were less biased than the original observations and are merely a subsample of a larger set. The resulting rattlesnake observation data comprised 675 observations.

Environmental Covariates
We examined candidate and proxy environmental variables that best explained western rattlesnake distribution based on rattlesnake ecology and factors known to influence their distribution and abundance [35,39,56,75]. Except for climate variables, all of the environmental features were in ESRI raster format and had the same geographic extent and pixel size (30 × 30 m- Table S1). Because rattlesnakes rely on behavioral thermoregulation to maintain their body temperature [79,80], solar radiation and weather variables are relevant for predicting their occurrence [36]. Therefore, we used minimum monthly temperatures and mean monthly precipitation from the active (non-hibernation) period for rattlesnakes in BC (1 April to 30 September, 1 km pixel size) from Chelsa bioclim data available at: http://chelsa-climate.org/bioclim/ accessed October 2020 [81,82]. Note that although rattlesnakes are active in Idaho [68,70] and northern/coastal California and western Oregon [72] until October or November, this typically is not the case in BC. We created a raster layer to average the climate (temperature and precipitation) information from April to September. We then resampled this raster layer from a 1 km pixel size to a 30 × 30 m pixel size using the nearest-neighbor resampling method to keep all environmental data at the same resolution. We used long-term averages from the period of data available . Lastly, we built the 'Potential Annual Incoming Insolation' layer using SAGA GIS [83].
Ruggedness can not only impact snake mobility, but it can also create a diversity of microclimates enabling snakes to locate suitable thermal surfaces, even though the overall averaged thermal index may be unsuitable for a particular area [36]. Thus, a second key factor in estimating the distribution and abundance of rattlesnakes is topography [84,85].
Since many topographic features are highly inter-correlated (e.g., slope, topographic position) and may be used in combination as habitat cues by rattlesnakes, we used a 'ruggedness' index [36]. This was calculated directly using the 'Terrain Ruggedness Index' in SAGA GIS [83]. We also included elevation in models (upper limit 2000 m) from provincial Digital Elevation Models [86].
In BC, the strongest cover associations that rattlesnakes have are with rock outcrops, talus slopes, shrub-steppe/grassland, riparian and open pine forests [35,39,87]. Therefore, the third GIS layer that we used contained categorical land cover variables derived from a supervised satellite imagery Landsat classification from the North American Environmental Atlas [88] and updated provincial land cover [89]. We used all 19 categories of available land cover from the atlas and the province (Table S1) to make sure that no gaps in coverage existed and that the vintage of this land cover classification was the period 2009-2011. We also used a grassland land cover layer from the GIS sources of the Grassland Conservation Council of British Columbia (Table S1).
Roadways fragment snake habitat and are also a source of mortality [64,90,91]. The warm surfaces of roads are presumably attractive to snakes for thermoregulation and so many snakes are killed by vehicles, often intentionally (N. Bertram and J. Surgenor, personal communication and [92]). To test for the effect of transportation infrastructure as a predictor, we used a road layer for the province derived from GeoBC's Digital Road Atlas [93]. We included a range of road types, noting that road placement probably has a greater impact on rattlesnakes than road type. Road density was used as an index of the influence of roads on rattlesnake distribution and abundance and also may serve as a proxy for land cover fragmentation and quality of habitat for rattlesnakes. Similarly, other types of human disturbance could influence the distribution and abundance of rattlesnakes, including the low elevation parts of ski hills and golf courses which we derived from an agricultural land use database compiled by the BC Department of Agriculture (see Table S1).
While pesticide exposure (through bio-accumulation in prey items) and interactions with humans or machines can pose risks for snakes [65], agricultural landscapes also often harbor food resources, access to water and relatively stable vegetation types (e.g., orchards or pasture). We therefore included several agricultural land cover types in our models-including intensively-farmed cropland (comprised of cultivated crops, fallow land, transitional cropland), vineyards and berries (combined) and perennial agricultural land of ecological value known to be used by rattlesnakes (orchards-fruit trees; see Table  S1).

Model Development
ENMs are based on correlations between species' locations and environmental variables within the environment in which they occur [94]. We deployed ENMs to develop relative index of occurrence (RIO) maps that were used to predict occurrence of western rattlesnakes over BC. Based on their high performance (ROC area under the curve metric-AUC [95]), the ones we chose were Random Forests [96], MaxEnt [97], and Generalized Boosting Models, also called Boosted Regression Trees [98]; further details on why these model platforms were selected are in Supplement S1. Three ENMs were used in conjunction to reduce the effect of individual-model variation in model outputs (see [99]). We ran all models in the BIOMOD2 statistical package in R [100]. We also considered the recommendations of Hao et al. [101] and followed logical protocols (e.g., [102]). We finetuned the MaxEnt model parameters using the ENMeval package [103] and Random Forests and Generalized Boosting Models with the BIOMOD2 package (ntree and ntry parameters) [100]. We focus on the general patterns and areas where snake density was significantly higher ('hot spots') and areas where density was significantly lower ('cold spots') of the predictions.
Because many presence observations of rattlesnakes were opportunistic, we made a sampling bias layer (i.e., using the density of observation points). We then randomly selected 5000 pixels that had higher values in the sampling bias layer; thus, there was a higher probability of the pixel being selected as a background point. Using these background points in the presence/background model will result in them having the same bias as presence-observation points [104,105].
Having identified relevant environmental covariates as predictors for rattlesnake occurrence, we conducted exploratory analyses, including calculating Pearson correlations for all pairs of variables (Table S2). Removing highly correlated variables is not necessary for Random Forests, MaxEnt, or Generalized Boosting Models but we did this to avoid data redundancy, to ease interpretation of model outputs (i.e., response curves) and to speed up model computation time. Most variable pairs showed low correlations (< 0.50; Table S2); however, because there was a highly significant correlation between ruggedness and other topographic variables (e.g., slope, r = 0.89), we assumed that using only one variable would be adequate for the model and that relationships would be similar for both ruggedness and other topographic variables. Similarly, temperature and elevation were highly correlated (r = −0.95), and solar insolation and aspect were strongly related (correlations were not performed because aspect was a categorical variable). We considered elevation (DEM) and solar insolation to be the more proximate functional factors for rattlesnake habitat use [56,106], and therefore we omitted temperature and aspect from our final models. We compared model performance between Random Forests, MaxEnt and Generalized Boosting Models using the AUC for the most reasonable predictions and inferences [107,108]. We calculated metrics of model performance as performed by LeClere et al. [109] such as AUC [95] as described in DeLong et al. [110] and also included other metrics (Cohen's kappa [111] and true skill statistics (TSS) [112]). AUC provides a measure of the model's discriminatory ability (threshold-independent measures of model performance); thus, an AUC value of 0.5 demonstrates that a model's discrimination is no greater than expected by chance, whereas a value of 1 is indicative of perfect discrimination; values of >0.8 provide excellent discrimination [95]. Cohen's kappa statistic provides a correction for the overall accuracy of model predictions by the accuracy expected to occur by chance (see [113]. Thus, kappa values range from −1 to +1, where +1 indicates perfect agreement and values < 1 indicate that performance is no different than would occur randomly [111]. Lastly, TSS, corrects for the dependency problems of kappa, but retains all of the advantages of the former method [112]. Recently, estimation of kappa and TSS from presence-only models has been criticized for two reasons. First the lack of true absences precludes the computation of an unbiased confusion matrix, and second, in the case of kappa, it is overly sensitive to presence [114,115]. Therefore, to verify our results, we conducted a second model run for the ensemble model and computed the Boyce Index. This index provided an independent measure of model performance and is largely insensitive to species prevalence and provides additional information on model quality [116]. To evaluate the predictive power of our models, we followed Huettmann and Gottschalk [117] and divided data into two groups: 1) a training dataset (75% of data) used to create the initial model; and 2) an independent test dataset (25% of data) which we used to test model quality [95] and to ensure accuracy and repeatability. This approach focuses on predictions first and then inferences from those predictions [107]. Although it may not be an issue with our data, we controlled for the possibility of over-fitting by using a regularization multiplier (RM) of 1 (as in the default version of MaxEnt) which we considered appropriate for our model. We used the ENMeval package in R [103] and checked RM values ranging from 0.5 to 10, in 0.5 increments. The combination of RM values (20 different values) and feature classes (FCs, 5 different classes) resulted in 100 different model combinations from which we selected the best model using Akaike's Information Criterion (AIC). Lastly, we created an ensemble (averaged) model [99,101,118] from all model platforms which we used to calculate our population estimates. This was carried out using BIOMOD2 and was performed by overlaying the lattices from the models and based on their TSS values. To convert the continuous predictions (from 0 to 1) to binary values, we developed thresholds at the mean RIO value as well as values of 0.5, 0.6, 0.7, 0.7, 0.75, 0.8 and 0.99 (Table S3).

Population Estimate for BC
We used population density data from three intensively-studied small study areas in BC (Figure 1; see Table 1). The study sites included (1) Osoyoos-an 11-year study (2002-2012) based on intensive capture-mark-recapture (CMR) of snakes [28]; (2) White Lake Basin area, approximately 34 km north-northwest of Osoyoos, where densities were estimated between 2015 and 2018 [63,119]; (3) Kalamalka Provincial Park and neighboring ranch near Vernon, approximately 140 km north of Osoyoos (2018-2020) [51]. Note that we only included the most recent estimates for adults in 2020 from the ranch and provincial park in Vernon combined (mean 65.9, 95% CL 58-4-73.4) as there has been a 10-33% population decline in this area since 1983. We therefore calculated the adult population estimate as a product of rattlesnake habitat availability (km 2 ) times the density estimates from each of the Osoyoos, White Lake Basin and Vernon study areas and then taking an overall average. We then determined the threshold at which rattlesnakes would be considered present [120]. To do this, we first chose a critical RIO threshold value from the ensemble model from Random Forests, MaxEnt and Generalized Boosting Models. Based on prior estimates from den counts [121], we investigated several candidate thresholds (Table S3) and chose an RIO threshold value of 0.8, which was approximately 75% of the maximum, from sightings (see also [108]). This value provided the most realistic population estimate, was consistent with estimates from counts of snakes at dens by J. Hobbs [121] and was also verified by expert opinion (J. Hobbs, pers. comm. provided a provincial estimate of an upper limit of ~8K western rattlesnakes in BC and M. Sarell, pers. comm. inspected the ROC maps for accuracy). We divided the land cover pixels into two groups ('suitable' and 'unsuitable' land cover) by using predicted pixel values based on the sampled sighting locations at the 0.8 threshold. Pixels with values greater than or equal to this threshold value were considered suitable land cover (i.e., 'habitat') for rattlesnakes in our model outputs. We then converted the pixels of suitable land cover into polygons. We calculated how many pixels were included in the RIO model in the additional range area (km 2 ) equal to or above that threshold value. We extracted all land cover polygons within or intersecting the buffer polygons generated in the first step and used these as the basis of the total area of rattlesnake habitat in BC. We then calculated the area of suitable land cover within the buffered pixels. We also used a very conservative threshold to calculate our final population estimate, which meant that only the land cover areas with the highest suitability were included (see thresholds; Table S3).

Scenarios of Protected Areas, Roads and Agricultural Land
We ran three scenarios to evaluate the impact of: (1) adding or subtracting protected areas, (2) adding or subtracting the area of roads (buffered-see later), and (3) adding or subtracting agricultural land on the western rattlesnake population estimate. For each scenario, the area of suitable land cover for western rattlesnakes was calculated (km 2 ).
In our first scenario, we counted pixels included in the RIO threshold within protected areas in the distributional range and those outside. The protected areas layer (Figure 1), which included National Parks, National Wildlife Areas, Migratory Bird Sanctuaries, and provincial parks, was obtained from Environment and Climate Change Canada [122]. While on the one hand protected areas could protect western rattlesnakes from persecution, on the other hand some protected areas have high numbers of human visitors and disturbance can influence movements or survival of snakes.
The second scenario we ran included and excluded roadsides from the population estimate calculations. We did this by buffering and masking a 500 m area from roads (Figure 1). The 500 m criterion was based on a study demonstrating that timber rattlesnake (Crotalus horridus) and Louisiana pine snake (Pituophis ruthveni) numbers were depleted by > 50% in areas occurring within 450 m of roads with moderate traffic [123]. Note that we were not able to determine traffic volume from GIS and therefore selected mostly paved roads, which we assumed, on average, would have moderate traffic and perhaps pose a greater threat of road mortality to snakes [124].
Depending on farming intensity and cover types, agricultural land may or may not be a suboptimal land cover type for rattlesnakes in BC [48]. Therefore, for our last scenario, we postulated that rattlesnakes avoided agricultural lands (Figure 1). To test the influence of agricultural land on the RIO models and subsequent population estimates, we calculated the number of pixels assigned to land cover suitable for rattlesnake occupancy inside and outside of agricultural land.

Relative Index of Occurrence (RIO) Models
In order of performance, Generalized Boosting Models were best (AUC 0.  (Table 1; we also compared estimates using each study region separately). The separate ensemble model run to cross check results using the Boyce Index indicated similarly high performance: ensemble model, 0.95; GBM, 0.90; Random Forests, 0.84, and MaxEnt, 0.90.
The RIO maps indicated that the highest-quality land cover occurs in many of the same areas as known dens for rattlesnakes, providing confirmation of the accuracy of our models. In the Thompson-Nicola range (Figure 1), the highest-quality indices were in the central portions, while in the Okanagan they are often along edges (lakes, such as Okanagan Lake) and in Vernon high concentrations were at the western edges and southern part of the range (Figures 2 and 3). Examination of the mean values for each of the small study areas generally demonstrated that Osoyoos had the highest RIO values, followed by Vernon Provincial Park, Vernon Ranch and White Lake Basin (Table 2). For all models, the RIO values were much lower for the entire range, than within the small study areas ( Table 2). This was expected as the small study areas had high habitat quality for western rattlesnakes.  (Figure 1a,b). The closer the RIO values are to 1, the higher the probability of occurrence of the species. Note the scale is in meters.  (Figure 1a,b). Note the scale is in meters.
The most important variables in all model platforms were elevation followed by annual insolation, land cover and ruggedness (Table 3). When the three model platforms were compared, there were some differences in the percent contribution of the variables ( Table 3). Inspection of model response (partial dependence) curves for relationships between rattlesnake presence and different environmental parameters matched existing biological knowledge for the species ( Figure S1). Overall, western rattlesnakes demonstrated a positive response to elevations between 0 and 1000 m, and increased with higher solar insolation and increased ruggedness (to a plateau).

Range-Wide Rattlesnake Population Estimate in Canada and Scenarios
Although the range of estimates depended on the RIO threshold value used (Table  S3), the final one chosen based on prior den counts was 0.8, giving a very realistic mean estimate of 9722 (±SD 3009) adult snakes from the ensemble model (separate population estimates using local densities from Osoyoos, White Lake and Vernon provided estimates of 10,702, 6344 and 12,118, respectively; see Table S3). The overall population estimate using the mean RIO threshold of 0.446 was extremely high (Table 4) and over-represented the suitability of land cover throughout the species' range (Table S3). We estimated the area of land represented by protected areas within the western rattlesnake's range in BC to be 923 km 2 . Within these protected areas, only 21.63 km 2 (2.3%) was available as the most highly suitable land cover, giving an estimate of 1144 (±354) adult snakes (or 11.8% of the total population estimated) within protected areas. A large proportion of the most highly suitable land for rattlesnakes was within 500 m of roads and estimated at 170.6 km 2 . Taking the mean estimate from densities in the three small study areas therefore gave a total of 9017 (±2791) adults (92.8% of the total population). Although agricultural land area within the rattlesnake range in BC was estimated to cover an area of 298 km 2 , the area of farmland identified by our models as being of high suitability for rattlesnakes was extremely small (<1 km 2 ).

Discussion
Our results provide the first quantitative and spatially explicit population estimates for western rattlesnakes across its Canadian geographical range using observational data and capture-mark-recapture results. This study of model-based estimates differs from some previous studies in that it was based on the range-wide data from BC (Conservation Data Centre), GIS land cover, topography, land use and climate data, much of which are publicly available (albeit with restrictions in place for access to rattlesnake observational data due to the threatened status of the species). Other recent studies of edge of range rattlesnake populations using modelling approaches include those by Hecker et al. [72] who used MaxEnt to develop ENMs in California and Oregon and Yousefi et al. [126] who conducted similar research in Iran. However, our approach is novel in that it combines results from intensive studies, including an 11-year demographic study at the Osoyoos Indian Reserve [28] and shorter-term studies at White Lake [63,127] and Vernon [51], with observational data of western rattlesnakes from throughout its range in BC. Our models reveal that only a small proportion of the population is likely conserved within protected areas. Our models also provide support for concern regarding the potential direct and indirect impacts that roads have on rattlesnakes in BC due to the large area of suitable land cover adjacent to major roadways in valley bottoms.
Our main goal of confirming previous den count estimates was supported by the estimated 9722 adults from our modeling and is probably the most conservative, realistic and defensible. One of the criteria for federal designation of the western rattlesnake in Canada as 'threatened' requires an estimated population size of <10,000 individuals-this specific criterion was not invoked in the most recent re-assessment of Western Rattlesnakes [46], but rather the numbers and decline of mature individuals were cited. An earlier approximation of the number of western rattlesnakes in BC suggested that there "could be fewer than 2000 or over 25,000" total individuals (COSEWIC, 2004 [46]). This estimate was based on Macartney's (1981Macartney's ( -1983) den counts (8 to 133 individuals per den) at a specific study site in the north Okanagan, and a provincial estimate of 392 confirmed dens. More recently, Hobbs [121] used head counts at 318 dens to arrive at an estimate of 3943 to 7896 rattlesnakes in BC. Despite the fact that estimating the total population from single or few den visits is subject to detectability biases-not all snakes present at dens can be counted and many dens are hard to find-the upper limit of the estimate from dens is aligned (~20%) of our estimates and given that both methods have inherent error. Therefore, we believe that our population estimates for adults based on quantitative mark-recapture estimates and most of the available rattlesnake observational data confirm the most robust estimates available to date. The estimate from the mean threshold may represent a potential population size for the western rattlesnake in BC, based on land cover and other environmental covariates.
Our overall adult population estimates could be inflated due to the fact that we used density estimates from the Okanagan valley, which has been identified as having the land cover types most highly correlated with rattlesnake presence in Canada and has the longest active season [46]. We drew upon data generated from three intensive rattlesnake studies, all occurring in the Okanagan, thus potential inflation of our estimates could have occurred given that we extrapolated to the rest of the snake's range. For example, our density estimates from Osoyoos, where rattlesnakes have been monitored and conserved over the long-term [28] were high, despite the fact that the southern portion is experiencing extremely high rates of human development; in contrast, in the north Okanagan (Vernon), suitable rattlesnake land cover is relatively undisturbed. The other intensive study area, White Lake Basin, is protected and essentially 'pristine' [63], harboring relatively high densities of rattlesnakes that are still considerably lower than that at Osoyoos (34.5/km 2 ). Both densities at the Osoyoos and White Lake are much lower than those at Vernon (see range of estimates in Table S3). Although, Macartney [87] found overall densities of 155.8 individuals/km 2 in the two neighboring subpopulations (provincial park 123.4/km 2 ; ranch 219.6/km 2 ) near Vernon, many changes have occurred in this area over the intervening 30 years. Some of the dens in the Vernon study are no longer occupied, and the population has declined by ≈50% (overall density 93.5/km 2 , provincial park 61.9, ranch 151.7 [51]). Although taking a mean estimate from all three study areas is conservative, nevertheless within the BC range of western rattlesnakes these three areas have high habitat quality.
It is also important to note that our models did not simply extrapolate densities from localized estimates of snake densities to the entire range within BC. The purpose of the RIO map was to control for land cover suitability in the rest of the range and we purposefully included only cells with high 'habitat quality' in our population estimates. As with all modelling exercises, our outputs are subject to potential biases due to the limitations of the data that we incorporated. For example, the only climate data available were at 1 km resolution and we had to rescale these data to 30 m resolution to match our other predictor variables. We do not think this will have a great influence on our model results for three reasons: (1) climate variables would not change significantly at finer scales but are important to include; (2) including topographic variables such as ruggedness at a fine resolution is important because of considerable topographic variation within the areas used by rattlesnakes; using a 1 km resolution would lose much of this critical detail; and (3) the rattlesnake observations were a mix of more recent GPS-based point data (3 m accuracy, with a buffer of 10-12 m), or derived from the centroids of larger variable-sized polygon locations (almost all 250 m accuracy).
All models and estimates of population size involve assumptions, and require decisions to be made, but their products are more defensible and repeatable than estimates based on anecdotal observations or expert opinion alone. Because the Conservation Data Centre and other observations of rattlesnakes were not randomly sampled, spatial survey biases could occur by over-representing highly sampled areas (e.g., provincial parks or other protected areas) or because the data were not statistically independent (spatial autocorrelation due to multiple observations from the same individuals).
Although our ENMs models performed well, and high AUC values were obtained, it should be emphasized that these models represent a first approximation. They are also subject to sampling bias since incidental observations of rattlesnake housed by the Conservation Data Centre are more likely to be gathered from accessible areas such as those close to roads or residential areas. Supporting this, a large percentage of the highest-quality land cover from the models, and consequently a large proportion of the estimated population, was found within 500 m of roads. Whether this is a real effect (roads pass through the valley bottoms which represent the highest-quality land cover) or an artefact due to sampling biases (more snakes are observed close to roads in the opportunistic Conservation Data Centre observations) requires further research. However, it is known that almost all habitat for rattlesnakes in BC is bisected by roads so controlling for this sampling bias is extremely difficult [121]. Note that we did attempt to minimize potential biases in the background data by building a sampling bias layer (see methods and details in [104,105]. One area for consideration is further modification (ideally using LiDAR data) by stratifying within suitable vegetation types and ranking specific crown closures and vegetation classes.
It may be useful to further separate the data into different subsets and build specific ENMs for training and testing within these subsets. For example, separate models could be built for dens or subsets of older versus more contemporary observations. Separate models for historical (e.g., pre-2002 or pre-GPS availability in 1997) versus contemporary (post-2002, aligning with individual study site temporal spans) may provide a handle on the influence of population declines on our models. Although we did conduct separate models for dens in earlier model runs, we could only use a limited subset of the approximately 477 dens (392 confirmed) in the province because of the confidentiality reasons previously mentioned (see Materials and Methods). To confirm our results, future models could consider using abundance data from these den counts in relation to the landscape mosaic and connectivity for rattlesnake habitat using point process models (e.g., [128]). Using dens as observation data points for rattlesnakes must also consider that reported den locations may be biased as not all dens are found.
To improve future models, more detailed analysis using demographic models, as well as their assumptions, effects of landscape heterogeneity, climate data and solar radiation where data become available would be useful. We used spatially explicit estimates of presence based on modelling which were adjusted according to the intensive study area average estimates. More process-models such as spatially explicit population viability models could also be developed [129].
Our estimates suggest that although the species is not imminently threatened with extirpation because there are ~10,000 estimated individuals, there is cause for concern about the long-term viability of the population [127]. This is especially so given effects of habitat loss and degradation and the likelihood that, in its Canadian range, the western rattlesnake population is declining [28]. This situation, coupled with the fact that rattlesnake habitat loss continues while incorporating regularly updated satellite land cover imagery presents challenges and makes estimates of this species' population size a 'moving target'. For example, Atkins [51] recently demonstrated that the population at Vernon has declined by 10-33% over the last 30 years. These population declines in the western rattlesnake may parallel those in another large snake in BC-the Great Basin gophersnake, (Pituophis catenifer deserticola) [130]. For that species, models estimated adult population declines of up to 40-50% over three generations based on road mortality and other threats such as accidental and intentional killing [131]. Similarly, western yellow-bellied racers (Coluber constrictor mormon) experience high road mortality, even in comparison to rattlesnakes and gophersnakes [63].
Given the status and apparent decline in abundance of western rattlesnakes, it is of concern that (according to our models) a low percentage of the population (based on highest-quality land cover) is projected to occur within protected areas. A similar situation occurs in the case of the Great basin gophersnake, which also has its core geographic range in Canada in the Okanagan valley [132]. Supporting this, Deguise and Kerr [133] found that only 1.5% of areas with > 30 endangered species were protected in Canada. It is therefore critically important that non-protected lands be managed for rattlesnake conservation. For example, extensive rangelands appear to provide habitat for rattlesnakes and in some areas (e.g., Vernon) have higher densities than a neighboring protected area (see above [51]). There is also evidence that rattlesnakes in provincial parks, or other areas with high human use, rattle less frequently due to disturbance [51]. Our results support findings elsewhere that protected areas alone are not sufficient for biodiversity conservation, and as in the case of the western rattlesnake, their value is dependent on representation, access and human visitation levels, and management effectiveness [134].
The most contiguous area of suitable land cover for the western rattlesnake in BC occurs in the most southerly portion of its range where agricultural land and roads are most prevalent (i.e., the Okanagan valley). Over the past 100 years, 77% of Idaho Fescue (Festuca idahoensis spp. idahoensis) and 68% of the antelope-brush (Purshia tridentata) and needle-and-thread grass (Hesperostipa comata)-important vegetation types for western rattlesnakes-have been lost and replaced by agriculture and urban land use in the Okanagan [135,136]. Geographically, conversion to agriculture or urban and industrial uses has resulted in more than 47.6% loss of native grasslands in the northern Okanagan Basin and substantial losses in the South Okanagan Highlands (38.6%), South Okanagan Basin (20.5%) and Thompson Basin (20.0% [52]). Recent trends for land conversions to vineyards in the southern interior of BC has further reduced land cover quality for rattlesnakes by removing previously heterogeneous agricultural landscapes containing hay fields, orchards and vegetable crops, interspersed with natural vegetation. As well, preparation for vineyards includes removal of topsoil and subsoil contouring, removing the natural cover that harbors rodent prey populations. At least some agricultural areas could be attractive to snakes because they provide water, diverse landscapes as well as prey, and contain patches of natural or semi-natural land cover [7]. Although Spear et al. [56] did not include the area of agricultural land in models directly, they determined 'resistance' values for land cover types, including agricultural land, which were assigned the same value as invasive vegetation, railways, low housing density or very rough terrain. Recent evidence suggests that rattlesnakes in disturbed landscapes have lower body condition than those in undisturbed landscapes in BC [75]. Thus, agricultural land and urban areas may act as population sinks, as may areas where road expansion is ongoing [123]. If managed in a holistic way (such as organic farming regimes, including intercropping with native vegetation), agricultural land may provide habitat for rattlesnakes and snakes should be conserved for the ecosystem services they provide (e.g., preying on agricultural pests such as mice and voles). Clearly, more research is needed to investigate the role of agricultural land in supporting or impacting this species, especially given the rate of land conversion. Because of the coincidence of rattlesnake distribution with areas under residential development in BC and/or due to their attraction to roads, our study highlights and further supports the potential negative impact of roads because of the close proximity of high value rattlesnake land cover types. Indeed, road placement may have a greater bearing on road mortality than road type, given that many roads are in valley bottoms where most snakes occur [121]. Even at low traffic volumes, roads in BC can have a devastating impact on rattlesnake populations, as demonstrated by the recent Population Viability Analysis by Winton et al. [127].

Conclusions
Estimating population sizes of species across their range (within and across political jurisdictions) is important for setting targets and evaluating population changes over time. This can be used as a baseline for ongoing monitoring, management and species conservation. Similar studies to ours have been carried out for other rare species-scalysided mergansers (Mergus squamatus) and various crane species in China [137,138], as well red pandas (Ailurus fulgens) in Nepal [108]. Our models serve as only one important approach for trend monitoring and future re-assessments of the status of western rattlesnakes in BC. They may not be sufficiently sensitive to monitor trends in a declining population. However, models based on new observations of rattlesnakes, ideally including density estimates from a broader area, could be used to verify and track potential declines, especially if paired with data on population structure and recruitment. Not only do models such as ours provide critical information for the conservation and management of the western rattlesnake in BC but they also can be used for spatial conservation planning purposes and wider conservation decision making [12,139,140]. Most importantly, our models provide the first quantitative and spatially explicit population estimates for western rattlesnakes across its Canadian geographical range.
Models such as ours can be used to quantify the impact of habitat loss and other threats for landscape-level conservation planning purposes and habitat restoration using Marxan (https://marxansolutions.org/ accessed 21 September 2021), prioritzR [141] or other spatial prioritization tools and/or cumulative effects assessment [108,[142][143][144][145][146][147][148]. Our study provides an empirical framework, including a transparent and repeatable process, to conduct such models for other species at risk and their habitat. If these were combined and overlaid, then the identification of hotspots for species at risk would allow threats and protected areas (e.g., critical habitat for recovery) to be prioritized spatially and in relation to competing resource development, thereby facilitating persistence of a broader group of at-risk species' populations.
Supplementary Materials: The following are available online at www.mdpi.com/article/10.3390/d13100467/s1, Table S1. Descriptions and sources for variables included in models. Supplement S1. Model selection. Table S2 Supplementary Results: Pearson correlation coefficients for environmental covariates. Table S3. Thresholds. Figure S1. Response curves of rattlesnake occurrence in response to predictor variables used in individual models and ensemble model. Each graph shows how rattlesnakes' Relative Index of Occurrence changes based on different values of each variable.  Institutional Review Board Statement: Not applicable (no animals were handled during this specific study).

Data Availability Statement:
The rattlesnake observation data are archived by the British Columbia Conservation Data Centre as part of the Wildlife Species Inventory (WHI) managed by the BC Ministry of Environment. Because the species is classified as Threatened under the Species at Risk Act in Canada (2003), and suffers from direct and indirect persecution, these data are under strict confidentiality agreements; contact Conservation Data Centredata@gov.bc.ca. Digital elevation model spatial data are available upon request from the BC Ministry of Forests, Land, Natural Resource Operations and Rural Development. The Agricultural Land Cover data are available on special request from the BC Department of Agriculture (https://catalogue.data.gov.bc.ca/dataset/92e17599-ac8a-47c8-877c-107768cb373c-accessed 4 November 2020. Land cover data are freely available at http://www.cec.org/files/atlas/?z=3&x=-93.1641&y=61.9803&lang=en&lay-ers=polbounds%2Clandcover2015ls&opacities=100%2C100&labels=false -accessed 10 October 2020; climate data at https://chelsa-climate.org/-accessed 10 October 2020; Grassland Conservation Council of British Columbia GIS grassland land cover; https://a100.gov.bc.ca/pub/acat/public/viewReport.do?reportId=54174; accessed 23 November 2020; DEM data from the Government of British Columbia at https://www2.gov.bc.ca/gov/content/data/geographic-data-services/topographic-data/elevation/digital-elevation-model-accessed 4 November 2020; and road data from: https://www2.gov.bc.ca/gov/content/data/geographic-data-services/topographic-data/roads-accessed 10 October 2020.