Factors Affecting the Detection of an Imperiled and Cryptic Species

Population surveying and monitoring are important for identifying conservation needs and tracking trends in populations, communities, and ecosystems over time and laying the groundwork for conservation management and policy decisions. If species or populations go undetected because of inadequate effort or sampling design, protection and management cannot be properly provided. Due to the widespread loss of populations, the Eastern Massasauga (a rattlesnake) was recently listed as a federally threatened species in the United States; it is also listed as threatened in Canada. Given its current conservation status, there is considerable interest at state and federal levels in determining how to best survey for Eastern Massasaugas to aid in management decisions. Using a 16-year dataset, we examined the relationships among environmental, temporal, area, management, and search effort factors on the detection probability of Eastern Massasaugas. We found that four abiotic parameters (solar irradiance, shaded air temperature, three-day maximum air temperature, and humidity) and three search parameters (effort per researcher, search area, and search time of day) influenced detection of Eastern Massasaugas. As the current biodiversity crisis continues, the cost-effective use of resources and scientific expertise will continue to increase in importance. We hope our results stimulate similar analyses in other taxa, which will be critical for designing and implementing regional survey and monitoring programs.


Introduction
Surveying and monitoring are important for tracking temporal trends in populations, communities, and ecosystems, and they lay the groundwork for conservation management and policy decisions [1]. However, surveying for rare, imperiled, and/or cryptic species can be a challenging task [2,3] because a large amount of effort may be required to determine presence/absence at a single location, let alone across a species' range [4]. If inadequate sampling design or effort result in false negatives, populations may not receive the protection and management needed to persist [4,5]. Consequently, determining the appropriate design and effort required to evaluate status and trends of imperiled and cryptic species is a vital first step for surveying and monitoring programs [6,7].
Reptiles are undergoing global declines due to habitat loss and alteration, disease, global climate change, environmental pollution, and overharvesting for the pet trade [8][9][10]. Over the past decade, the number of reptiles on the IUCN (International Union for Conservation of Nature) Red List doubled,

Study Species and Sites
The Eastern Massasauga (Sistrurus catenatus) is a small rattlesnake ranging from southern Ontario and central New York west to Iowa and Missouri [29]. The geographic distribution of extant localities is fragmented because of range-wide extirpations [30]. Thus, the Eastern Massasauga is listed as endangered or threatened in every state it occurs except Michigan, where it is a species of special concern [30]. Recently, the Eastern Massasauga was listed as threatened under the United States Endangered Species Act of 1973 [31,32]. In Canada, they are listed as threatened under Canada's Federal Species at Risk Act [33]. Given its current conservation status, there is considerable interest at state and federal levels in determining how to best survey for Eastern Massasaugas to aid management decisions.
We conducted our study at Carlyle Lake (Clinton County) located in south-central Illinois, USA (38.62 • north (N), 89.35 • west (W)). Carlyle Lake is a~10,500-ha human-made impoundment (1967). The lake is surrounded by~4500 ha consisting of remnant grasslands, wetlands, and upland and bottomland forests [34]. We conducted our searches in 13  Although S. catenatus do not communally hibernate like many other rattlesnakes, they do aggregate in wet grassland habitats and over-winter in crayfish burrows. This artificial inflation of densities affords easier access to large numbers of individuals. Thus, all visual encounter surveys occurred during the emergence period which ranged from late February until mid-May depending on general climatic conditions (earlier spring weather resulted in searches conducted in late February, and late or cooler spring conditions extended searching into mid-late May). Searches consisted of surveyors spread approximately 2-3 m apart walking parallel transect lines through suitable grassland habitat. We constrained search areas between physical habitat boundaries such as roadways, woodland edges, and other urban and recreational land features. To vary the amount of effort spent per site, we attempted to vary search times by either increasing or decreasing the number of parallel transects run. During each survey, searchers would scan the ground, search in brush piles (as encountered), and intermittently move vegetation with snake hooks in areas where vegetation was dense. Upon encountering a snake, all searching and time recording stopped until snake-specific data was collected.

Environmental Variables
We used surveys conducted from 1999-2014 having complete data for the analysis. At the start and end of each search, we recorded shaded air temperature ( • C), substrate temperature ( • C), relative humidity (%), and wind speed (m/s) using a Kestrel®pocket weather meter. We then generated two additional variables representing the change in shaded air temperature and change in substrate temperature during the search period. From the local weather station at the US Army Corps of Engineers office (38.62 • N, 89.38 • W), we collected the previous day minimum, maximum, and mean air temperatures ( • C), as well as the previous three-day average of the minimum, maximum, and mean air temperatures ( • C). We calculated the area (ha) of suitable habitat searched using geographic information system (GIS) data layers, ground-truthing, and then summing the area searched. For solar irradiance, we used the modeled solar radiation on a horizontal surface received from the sky from the National Solar Radiation Database (http://nsrdb.nrel.gov). The irradiance data were derived from the Physical Solar Model ver. 3 at 4 × 4 km resolution around the approximate coordinates between all study sites (38.61-89.34). Finally, for use in the analysis, we aggregated the solar data to represent the maximum solar radiation by day. We also collected historical burn data for all sites from the land managers and coded the variable as burn (if a full or partially controlled burn occurred at the search site in a given year) and no burn (if no controlled burns occurred at the search site). For search effort related variables, we recorded the number of searchers and the total search time (minutes), and we then calculated the mean search time per searcher (min). Finally, we converted the date of searches into day of the year. Because all variables were on different scales, we standardized them using a z-transformation; thus, they are represented as standard deviates. Summaries and explanations of all variables used are presented in Table 1.

Multicollinearity Assessment
To minimize multicollinearity and its negative effects on parameter estimates, we firstly used pair-wise correlation analysis to view correlations among variables. Then, we used ridge regression to determine which parameter estimates were unstable given some perturbation (penalization) using the lm.ridge function in the MASS package of R [35,36]. Briefly, ridge regression imparts a penalty term (λ) on β parameter estimates. In situations with multicollinearity, if β parameter estimates show large changes (numerically or directionally), they are unstable. For example, if we observed multicollinearity in the suite of variables describing the previous day's temperature, we would only use one representative variable from the suite per model. Once we removed all offending variables, we performed a second ridge regression to ensure no unstable variables remained. If multicollinearity effects persisted, we iteratively repeated the process until our reduced variable had stable β parameter estimates. Our correlation analysis revealed substantial associations between the three-day and previous day mean, min, and max temperatures and warranted caution using them in the same candidate models ( Figure S1, Supplementary Materials). Our initial ridge regression analysis revealed the following parameters were unstable: eSAT, eSUB, dSAT, dSUB, eWind, and eHUM. We eliminated those variables from subsequent analyses. We resolved all parameter instability from our variable set in our second (and final) ridge regression analysis, and all parameter estimates remained stable regardless of the penalty used (Table S1, Figure S2, Supplementary Materials).

Statistical Analyses
Because Eastern Massasaugas occupied all our search sites, we could focus on factors affecting detection probability. Thus, we employed mixed-effects, binary logistic models using the glmer function in the lme4 R package [37] on the full dataset. Our dependent variable was the presence or absence of Eastern Massasaugas during the search, and our random effect was the site. Using our reduced variable set (multicollinearity analysis above), we constructed 25 a priori models including a global and null intercept-only model ( Table 2). We also included quadratic terms for some variables because of the potential to exhibit non-linearities (Table 2). We assessed model fit following an information theoretic approach [38] using the AICcmodavg package [39]. Because our goal was to present the best predictive model for detection probability, we also constructed one post hoc model based on the β parameter estimates with confidence intervals not bounding zero from the global model and one model for land managers using readily obtained variables. We used the effects package to examine the responses of detection probabilities to the covariates [40,41]. Finally, we used the piecewiseSEM package in R to calculate the conditional and marginal r 2 values for the two highest ranked models [42]. To assess the predictive capacity of our models, we estimated the area under the curve (AUC) and constructed response operator curves (ROC) for the top and global models for comparisons.

Conservation Applications
We used the Manager model to construct a tool in Microsoft Excel, which allows user input Burn, planned starting time (rTime), sSUB, m3Min, and mEffort. For flexibility, we allowed mEffort to be user input from nSearch and planned tEffort. We designed the detection tool to provide the estimates and confidence intervals for the detection/non-detection probabilities and their associated 95% confidence intervals. To obtain the number of searches, we modified the following equation from MacKenzie and Royle [28]: where p* is the of detecting a species at least once, p is the observed detection probability, and k is the number of surveys; we then solved for k as follows: We then applied the detection model to a historic Eastern Massasauga dataset to estimate the probability of a false negative (occupied site but no detections). Although the detection probability does not typically suffer from false positives with trained surveyors, true negatives (snakes are not present) and false negatives (snakes are present but missed) comprise the non-detection probability. For the analyses, we conservatively assumed that non-detection always represented false negatives (occupied sites) and, thus, the product of the non-detections becomes the overall probability of missing Eastern Massasaugas from every survey. The historical data we used comprised 58 search events over four years (2008,2009,2010, and 2015) for a grassland site in the Chicago, Illinois, USA, region (Table S3, Supplementary Materials). We generated m3Min data from the National Climate Data Center website (ncdc.noaa.gov) using the nearest weather station where we could get data for all dates (O'Hare International Airport, Chicago, Illinois, USA, 41.96, −87.93). We summarized solar radiation data from the same coordinates as the weather station data as above (Section 2.2). Finally, we took the Diversity 2020, 12, 177 6 of 17 product of the non-detection probabilities by year and all surveys to generate the probability Eastern Massasaugas were present but not detected in any of the surveys.

Data Composition
Our initial dataset comprised 895 individual search events from 1999-2014 across 13 hibernacula at four state/federally owned properties. We had 400 S. catenatus detections (44.6%), which illustrated that, although a rare species, our data neither had too few nor too many detections. We removed 182 surveys which were missing some covariate data. Many of these removed surveys were from the first few years of the study when we were iteratively improving our search design and, thus, added or removed variables. This narrowed the range of years from 2002-2014. The final data comprised 713 independent survey events (Table S2, Supplementary Materials) with most (632 surveys) occurring at the two IDNR properties (EHSP and SSSP) followed by the two USCOE properties (DERA and DWRA). We always placed maximal effort on SSSP because it also serves as our reference population for determining demographic patterns. Our final dataset had 315 S. catenatus detections across the 713 search events (44.2%), illustrating that the removal of survey events that were missing data did not impact the proportion of detections per survey.

Factors Affecting Detection Probabilities
Of the 27 mixed-effects binary logistic models we examined, the best model was the post hoc model, which carried an Akaike weight of 0.83 (Table 3). Variables included in the model were Burn, daRAD, mEffort, rTime, sSUB, m3Min, and all associated quadratic terms ( Table 4). The post hoc model explained a high amount of variance in detection probabilities (r 2 cond = 0.40; Table 3) with most of the variance explained by the fixed effects (r 2 marg = 0.28; Table 3). The Manager model excluding daRad also had strong support and carried an Akaike weight of 0.17 with a higher amount of variance explained (r 2 cond = 0.39; Table 3) mostly by the fixed effects (r 2 marg = 0.27; Table 3). The global model, with 30 parameters estimated, had little support and did not carry any appreciable Akaike weight; it was~15.0 ∆AIC C (AIC = Akaike information criterion) units away from the top model (Table 3). The global model explained only a bit more overall variance (r 2 cond = 0.47; Table 3), again, with most of the variance explained by the fixed effects (r 2 marg = 0.33; Table 3). For reference, the null model was the second to last model carrying no appreciable Akaike weight and being~102 ∆AIC C units away from the top model (Table 3). The post hoc model had a high AUC value (0.82, 95% confidence interval (CI) = 0.79, 0.85) and, thus, good predictive capacity, as well as the Manager model (0.82, 95% CI = 0.78, 0.84), whereas not much of that capacity was lost in the reduction of factors when compared to the AUC of the global model (0.82, 95% CI = 0.78, 0.84: Figure 1). The six variables comprising the post hoc model illustrates that we surveyed under a range of conditions (Table 4, Figure 2). Of our total surveys, 383 occurred in burned/semi-burned habitats and 330 occurred in burned habitats. Surveys were completed under conditions centered with a daRad ≈ 800 J/s·m 2 , rTime ≈ midday, sSUB ≈ 18 • C, m3Min ≈ 6 • C, and mEffort ≈ 50 min (Table 4; Figure 2).

Conservation Applications
Our detection probability tool offers all input variables from the Manager model (white cells) and allows the calculation of mEffort using nSearch and the planned tEffort for the event (Figure 4). To calculate the number of searches required, we allow the user to enter their desired probability of detecting an Eastern Massasauga (Figure 4). The tool will also produce warnings if a user attempts to enter values outside the bounds modeled ( Figure 4). All transformations and calculations are available, but not modifiable (gray cells), to the user on the second page of the workbook ( Figure S3, Supplementary Materials). In our example, we planned a search in burned habitat, starting at 1326 h, with a goal of 1 h of effort per searcher (4 total h). Given an sSUB of 16.2 and an m3Min of 12.5, we predicted our detection probability to be 0.67 (95% CI 0.40, 0.86). Thus, we would need 4-20 surveys under similar conditions to achieve a 95% chance of detecting an Eastern Massasauga in at least one survey ( Figure 4).
Initially, we had 58 surveys for the Chicago Region site; however, 12 had combinations of values with lower sSUB and daRad, and higher values of mEffort than we modeled and, thus, were removed from the analysis leaving 46 surveys (Table S3, Supplementary Materials). The overall probability presence but no detections for Eastern Massasaugas at the site varied from 2.56 × 10 −5 to 2.54 × 10 −28 by year (Table 6). When considering all surveys, the probability was 2.56 × 10 −46 (95% CI, 1.17 × 10 −83 , 3.90 × 10 −26 ) or, roughly, an average of a one in 1.64 × 10 45 chance (Table 6). Thus, if our model applies to other grassland habitats, and depending on the threshold decided by the governing land management agency, the site can likely be declared extirpated for management purposes.

Conservation Applications
Our detection probability tool offers all input variables from the Manager model (white cells) and allows the calculation of mEffort using nSearch and the planned tEffort for the event (Figure 4). To calculate the number of searches required, we allow the user to enter their desired probability of detecting an Eastern Massasauga (Figure 4). The tool will also produce warnings if a user attempts to enter values outside the bounds modeled ( Figure 4). All transformations and calculations are available, but not modifiable (gray cells), to the user on the second page of the workbook ( Figure S3, Supplementary Materials). In our example, we planned a search in burned habitat, starting at 1326 h, with a goal of 1 h of effort per searcher (4 total h). Given an sSUB of 16.2 and an m3Min of 12.5, we predicted our detection probability to be 0.67 (95% CI 0.40, 0.86). Thus, we would need 4-20 surveys under similar conditions to achieve a 95% chance of detecting an Eastern Massasauga in at least one survey ( Figure 4).
Initially, we had 58 surveys for the Chicago Region site; however, 12 had combinations of values with lower sSUB and daRad, and higher values of mEffort than we modeled and, thus, were removed from the analysis leaving 46 surveys (Table S3, Supplementary Materials). The overall probability presence but no detections for Eastern Massasaugas at the site varied from 2.56 × 10 −5 to 2.54 × 10 −28 by year (Table 6). When considering all surveys, the probability was 2.56 × 10 −46 (95% CI, 1.17 × 10 −83 , 3.90 × 10 −26 ) or, roughly, an average of a one in 1.64 × 10 45 chance (Table 6). Thus, if our model applies to other grassland habitats, and depending on the threshold decided by the governing land management agency, the site can likely be declared extirpated for management purposes.  Figure 4. Screenshot of the spreadsheet-based detection tool for grassland Eastern Massasauga surveys. Forecasted values are entered to estimate the detection/non-detection probabilities and the number of searches required to achieve a threshold level of detection assuming a site is occupied.

Discussion
Sit-and-wait foraging snakes (such as the Eastern Massasauga), which are inactive for long periods, might not be efficiently sampled using traps or artificial cover objects [43,44]; however, there is some indication that passive sampling can provide supplementary data [45]. Relatively inactive species are best sampled using time-and/or area-constrained searches. However, determining parameters influencing detectability becomes of paramount importance in the design and implementation of survey and monitoring programs. We found three abiotic parameters (solar irradiance, substrate temperature, mean three-day minimum air temperature), two search parameters (effort per searcher and time of day), and one management related parameter (burns) which significantly influenced detection of Eastern Massasaugas.

Modeled Effects
Detection probabilities were slightly higher for surveys in burned habitat (either full or partial). Burning is typically used as a management tool in Eastern Massasauga habitat to maintain early successional vegetation [46]. Eastern Massasaugas do not show alterations in movements, home ranges, or habitat use associated with burns, but do select habitat characteristics favoring thermoregulation [47]. Thus, they are likely more observable thermoregulating during the spring emergence period from greater distances in burned compared to unburned habitats. Timing or pairing visual encounter surveys (VES) with management burns can be an effective tool for increasing the chances of detection. However, caution must be exercised when burning after snakes emerge, which could cause significant mortality [46,47].
Both substrate temperature and mean three-day minimum air temperature influenced detection probability of Eastern Massasaugas in our study system. Such thermal relationships are expected for an ectothermic vertebrate. Harvey [48] found that Eastern Massasaugas in a Canadian population Forecasted values are entered to estimate the detection/non-detection probabilities and the number of searches required to achieve a threshold level of detection assuming a site is occupied.

Discussion
Sit-and-wait foraging snakes (such as the Eastern Massasauga), which are inactive for long periods, might not be efficiently sampled using traps or artificial cover objects [43,44]; however, there is some indication that passive sampling can provide supplementary data [45]. Relatively inactive species are best sampled using time-and/or area-constrained searches. However, determining parameters influencing detectability becomes of paramount importance in the design and implementation of survey and monitoring programs. We found three abiotic parameters (solar irradiance, substrate temperature, mean three-day minimum air temperature), two search parameters (effort per searcher and time of day), and one management related parameter (burns) which significantly influenced detection of Eastern Massasaugas.

Modeled Effects
Detection probabilities were slightly higher for surveys in burned habitat (either full or partial). Burning is typically used as a management tool in Eastern Massasauga habitat to maintain early successional vegetation [46]. Eastern Massasaugas do not show alterations in movements, home ranges, or habitat use associated with burns, but do select habitat characteristics favoring thermoregulation [47]. Thus, they are likely more observable thermoregulating during the spring emergence period from greater distances in burned compared to unburned habitats. Timing or pairing visual encounter surveys (VES) with management burns can be an effective tool for increasing the chances of detection. However, caution must be exercised when burning after snakes emerge, which could cause significant mortality [46,47].
Both substrate temperature and mean three-day minimum air temperature influenced detection probability of Eastern Massasaugas in our study system. Such thermal relationships are expected for an ectothermic vertebrate. Harvey [48] found that Eastern Massasaugas in a Canadian population were more detectable at intermediate body temperatures (20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) • C) and when they were more visible (>75% body exposure). Furthermore, he found that weather conditions which expressed thermoregulatory behavior (e.g., slightly cooler, sunny days) for intermediate body temperatures also resulted in increased detectability.
In contrast to other studies of snakes in Europe and North America, we did not find a significant positive effect of day of year on detection probability. In the Jura Mountains of Europe (France and Switzerland), Kéry [4] found a unimodal relationship between day of year and detection probability for three species of snakes (Asp Viper Vipera aspis; Smooth Snake Coronella austriaca; Grass Snake Natrix natrix), and the same unimodal relationship was found in North America for both Banded Watersnakes (Nerodia fasciata) and Black Swamp Snakes (Seminatrix pygaea; [49]). However, the Banded Watersnake and Black Swamp Snake were primarily sampled with passive trapping throughout a full season. Still, it is possible that such a trend may have emerged for Eastern Massasaugas if we sampled throughout the active season. Our data do suggest that fluctuations in environmental temperatures as the day of year progresses can result in successful search periods outside of the spring emergence period, and season and/or day of year might not always be a good proxy for environmental temperatures.
Detection probabilities of Eastern Massasaugas decreased as solar irradiance increased, and a similar pattern was also observed for the time of day. Although the two metrics are not significantly intercorrelated, solar irradiance peaks at astronomical noon, as well as throughout the activity season (strongest during the peak of summer), but it becomes uncoupled in overcast situations. During the spring and early summer, Eastern Massasaugas bask in the morning [29], which explains our increased detection rates early in the day. As temperatures increase, there is less need to bask; thus, they become more crepuscular and cryptic during the hottest periods of the year. Our decline in detection rates illustrates the reduced need for basking during periods when shaded air temperatures are greater than 24 • C. A similar trend was found by Erb et al. [7] for Eastern Box Turtles (Terrapene carolina), which had significantly higher detection rates earlier in the day. However, the pattern does not hold for all reptile species. Common Skink (Oligosoma polychroma) detections in New Zealand were the highest and least variable at ambient temperatures (12-18 • C) during daytime hours and did not decline until the evening [50].
We found that detection increased as search effort within a visit increased, but there were diminishing returns. Survey duration within a visit can significantly influence detection values in some, but not all, reptiles [5]. Furthermore, increasing survey duration was found to increase the power to detect raptor migration trends [51] and increase the number of species detected in anuran call surveys [52]. While search effort is incorporated into detection estimates, the effects of search area size are less commonly investigated. For Eastern Massasaugas, we did not find a relationship between search area size and detection probabilities. Our result is similar to Kéry [4], who found no effect of search area on the detection of three common snake species in Europe, although he did find that detection rates increased as population sizes increased. For relatively common species, one would expect higher detection rates when population sizes are larger since there is a greater chance of encountering the organism. Endangered species often have reduced numbers across the landscape and, therefore, can have small population numbers regardless of patch size. In this instance, if population size remains relatively constant across patches, but patch sizes increase, we would expect detection rates to decrease as the survey area increases [53]. Therefore, the next step might be to examine the interrelationship between density/abundance, patch size, and detection probabilities. Overall, our results highlight and reinforce the need to maximize search effort for endangered species.
Numerous other survey methods were used across the Eastern Massasauga's range, including road cruising, drift fences, and cover board arrays. Researchers were effective at encountering enough individuals using drift fence-funnel trap arrays for radio-telemetry studies of Eastern and Western Massasaugas [46,47]. The only study to directly compare multiple methods found that effort required (in person-hours) is lower for the more passive methods of carpet squares, cover boards, and drift fences [45]. The study also found sex-specific behavioral biases with females more likely to be encountered using coverboards and males more likely to be encountered with drift fences [45]. The overall results do recommend the use of multiple encounter methods for monitoring projects [45], and further work could be done to examine how detection probabilities differ among methods.
As MacKenzie et al. [24] pointed out, rare and imperiled species are simultaneously the species most in need of information on state variables (e.g., occupancy, abundance) and vital rates, and the species for which such information is the most difficult to obtain. Imperfect detection probability (i.e., the probability of detecting a species when present is less than one) must be accounted for when estimating state variables [24,54]. While detection rates can covary with predictor variables, detection probability often is treated as a nuisance variable in the occupancy or abundance modeling process [55,56]. However, by focusing on detection probability, researchers can evaluate the effects of sampling covariates such as abiotic factors, sampling method, and survey type [6,57,58]. Furthermore, this information can be used to determine the minimum number of visits required to a site to declare absence with a given level of confidence [4,27], which is a crucial first step in the development of survey protocols and monitoring programs.
The model we constructed provides good predictive capacity, but it was derived from the southern range limit in grassland habitat for a species with a historically broad distribution [29,59]. Over their broad range, Eastern Massasaugas occupy habitats ranging from wet to dry grasslands throughout the Midwest, to bogs, fens, and peatlands in the Great Lakes region, to open woodlands in northern Michigan and parts of Canada [29,59]. Therefore, not only may clinal variation exist but also variation due to habitat preferences. Thus, our model may only have the predictive capacity for grassland sites in the lower Midwest. We encourage the repetition of the methods and validation of the model across the species' range.
The only previous study examining detection in Eastern Massasaugas used 54 search events over two years from May-August paired with radio-equipped snakes [60]. From their 11 detections, they found detection probabilities approached 1.00 when the per searcher effort exceeded 90 min [60], similar to our study. However, they did not find the strong coupling with temperatures we found, and although they found that detection probabilities approached 0.80 on cooler mornings, confidence intervals were wide [60]. Future designs should also include stage classes as some stage classes may be more readily detectable than others as males tend to move more, gravid females tend to bask more often during gestation, and neonates/juveniles are extremely cryptic.

Conservation Applications
Giving a 0.10 spread from the maximal detection probability (in parentheses) up to a 0.95 probability for each variable, our data for Eastern Massasaugas indicate that, to maximize detection efficiency, surveys should be conducted preferably in burned habitat, between 9: To simplify the decision-making process for Eastern Massasauga surveys, we developed a spreadsheet-based detection tool as a guide. A researcher or land manager can enter forecasted values for the significant parameters that influenced Eastern Massasauga detection and determine the probability of detection for the upcoming survey, as well as the mean number of searches required. Additionally, one can estimate detection probability after the survey by entering the actual parameter values recorded during the sampling period. Such data could be used to determine if a site should be declared extirpated for management purposes. However, we do caution the use of the model outside of grassland habitats and, potentially, the region itself until it is further validated range-wide for the species. Given the decline of numerous snake species around the world, there is an urgent need to establish current status and distribution states for species of conservation concern. To deploy resources (i.e., time and money) in a cost-effective manner, and avoid wasted effort to establish presence or extirpation at a given site, we anticipate that our spreadsheet-based detection tool can be modified to meet the needs for other species.

Conclusions
As the current biodiversity crisis continues [61], the need for cost-effective use of resources and scientific expertise will be increasingly important [5]. When population declines are suspected, it is often difficult to separate true absence from inadequate sampling without incorporating detection probabilities; this is especially true for imperiled and/or cryptic species. However, in many cases, existing datasets can be leveraged to provide insight and starting points for imperiled species survey and monitoring programs. We encourage other researchers to use their datasets to examine factors influencing detection rates, which will be useful for designing and implementing regional survey and monitoring programs. These types of analyses will ensure that conservation and management agencies will have the necessary information to make informed decisions when surveying and monitoring for imperiled and cryptic species.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/5/177/s1: Table S1. Final results of ridge regression analysis showing the stability in parameter estimates for all variables retained in the detection analysis of Eastern Massasauga searches conducted at Carlyle Lake, Clinton County, Illinois, from 1999-2014; Table S2. Number of annual surveys conducted at each major site which had complete survey-specific data and were used in the detection analysis of Eastern Massasaugas at Carlyle Lake, Clinton County, Illinois; Table S3. Historical Eastern Massasauga search data for four years of searching at a site in the Chicago, Illinois, USA, region with the standardized search variables, predicted detection and non-detection probabilities, and associated 95% confidence intervals (Lower and Upper). All variable abbreviations can be found in Table 1; Figure S1. Pearson Correlation analysis of all candidate variables to be used for the detection probability analysis of Eastern Massasaugas from fieldwork at Carlyle Lake, Clinton County, Illinois during the 1999-2014 field seasons; Figure S2. Ridge regression plot to determine the multicollinearity of environmental and temporal variables used in analyses. Variable abbreviations can be found in Table 1; Figure S3. Example screenshot of the calculation sheet for the spreadsheet-based detection tool for Eastern Massasauga surveys. All calculated values are there for reference to see how calculations were made and are not editable; Excel workbook titled Final_Massasauga_Detection_Tool.xlsx as a companion guide for use in surveys.