Factors Affecting the Detection of an Imperiled and Cryptic Species

Abstract

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.

1. 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, and approximately 21% of reptile species are considered vulnerable to extinction [11]. Although many snake species are of conservation concern [12], conservation status assessments are hampered for over 80% of the species because they remain unevaluated or are data-deficient [11]. For example, in Illinois, United States of America (USA; where our study was conducted), approximately 40% of all snake species are currently listed as Species in Greatest Conservation Need [13].

We are only beginning to understand the varied roles snakes play in aquatic and terrestrial systems [14,15,16,17], but snake declines can have negative consequences on ecosystem structure and function [18]. Our general lack of information stems from their cryptic behavior, minimal activity patterns, and use of inaccessible (e.g., subterranean, arboreal) habitats [19,20,21]. Furthermore, such ecological traits result in a deficiency of high-resolution distribution and status data, which can obscure range-wide trends [22,23].

Establishing true presence/absence for non-imperiled snake species within a specified area can require substantial effort [4]. When population declines are suspected, it is challenging to separate absence from poor sampling without incorporating detection probabilities [24,25]. This challenge is particularly true for cryptic and imperiled taxa such as snakes, which led to increased interest in estimating and accounting for detection probability over the past decade [26,27]. With a known detection probability, one can determine the number of survey events required to detect a species with a given level of confidence [28]. Furthermore, incorporating survey-level covariates can help explain variation in detection estimates and provide insight into the factors that must be accounted for when surveying for a target species [27]. Once the number of surveys required to assess presence/absence is confidently established, researchers can then determine how many sites can be surveyed given the resources available (i.e., time, money). Prioritizing survey efforts during conditions associated with high detection can ensure accurate occupancy designation or improve abundance estimates. Such data are valuable for imperiled species.

Thus, the objective of our work studying a cryptic and imperiled snake was to determine the relationship among environmental, temporal, area, management, and search effort factors and the detection probability of the Eastern Massasauga (Sistrurus catenatus). Determining the factors associated with species detection aids in maximizing the amount of information per dollar spent on conservation efforts and provides an estimate of the number of searches required to achieve a certain confidence of true absence. Our results also provide further refinement of the optimal conditions and methods for conducting surveys to maximize survey efficiency, which will increase confidence in presence/absence determinations and reduce bias in the number of individuals encountered.

2. Materials and Methods

2.1. 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 grassland hibernacula patches at four public sites: (i) Dam East Recreation Area (DERA); (ii) Dam West Recreation Area (DWRA); (iii) Eldon Hazlet State Park (EHSP); and (iv) South Shore State Park (SSSP). DERA and DWRA are owned and maintained by the United States Army Corps of Engineers (USCOE), whereas EHSP and SSSP are owned and maintained by the Illinois Department of Natural Resources (IDNR).

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.

2.2. 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. Environmental, temporal, and effort variables used in the analysis of Eastern Massasauga searches conducted at Carlyle Lake, Clinton County, Illinois, from 1999–2014 with descriptions and units of measurement. Variables with a quadratic (squared) term included are represented by an asterisk.

Variable	Description	Unit
dSICA*	Presence or absence of an Eastern Massasauga in a search	binary
rTime*	Time calculated as the relative time of day on a 24-hour basis (hour + minute/60 + second/360)/24	h
DayofYear	Day of the year the search was conducted	count
sArea	The total area of grassland that habitat searches were conducted in/patch size	ha
Burn	Whether or not a burn was conducted in the search area	binary
nSearch*	Number of searchers per search	count
mEffort*	Mean effort per searcher (tEffort/nSearch)	min
tEffort*	Total search effort per search (cumulative time per searcher)	min
Precip	Precipitation for the day of searching	cm
daRad*	Historic solar irradiance data	J/s·m2
sHum*	Relative humidity at the ground level as the start of searching	%
eHum*	Relative humidity at the ground level as the end of searching	%
sWind	Wind speed ~1.5–2 m above ground level at the start of searching	m/s
eWind	Wind speed ~1.5–2 m above ground level at the end of searching	m/s
sSAT*	Shaded air temperature at the start of searching	°C
eSAT*	Shaded air temperature at the end of searching	°C
dSAT*	Change in shaded air temperature during searching (sSAT − eSAT)	°C
sSUB*	Substrate temperature ~2 cm in the substrate at the start of searching	°C
eSUB*	Substrate temperature at the end of searching	°C
dSUB*	Change in substrate temperature during searching (sSUB − eSUB)	°C
pMean*	Mean temperature for the previous day	°C
pMax*	Maximum temperature for the previous day	°C
pMin*	Minimum temperature for the previous day	°C
p3Mean*	Average mean temperature for the previous three days	°C
p3Max*	Average maximum temperature for the previous three days	°C
p3Min*	Average minimum temperature for the previous three days	°C

.

2.3. 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).

2.4. 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. Candidate model parameterizations used for the detection probability and capture yield analyses. All logistic detection models included a random effect of site.

Model	Variables	Description
Global	All Variables	Global additive
Null	Intercept	Intercept only
3DayTemps	m3Min, m3Max, m3Mean	Recent temperature
Burn	Burn	Habitat management
Calendar	rTime, DayofYear	Seasonality
EvapCooling	sHUM, sWind, sSAT, Precip, daRad	Evaporative cooling
m3Mean + daRad	m3Mean, daRad	Recent mean temps
m3Mean + Precip	m3Mean, Precip	Recent mean temps and precipitation
m3Min + daRad	m3Min, daRad	Recent min temps
m3Min + Precip	m3Min, Precip	Recent min temps and precipitation
tEffort + Burn	tEffort, Burn	Search effort and management
tEffort + sArea	tEffort, sArea	Search effort and area
tEffort + sArea + Burn	tEffort, sArea, Burn	Effort, area, and management
PDayTemps	pMin, pMax, pMean	Previous days temps
pMean + daRad	pMean, daRad	Previous mean and radiation
pMean + Precip	pMean, Precip	Previous mean temp and precipitation
pMin + daRad	pMin, daRad	Previous min and radiation
pMin + Precip	pMin, Precip	Previous min and precipitation
ProximateTemps	sSUB, sSAT	Proximate temps
SearchEffort	tEffort, nSearch, sArea	Total effort
sSUB + Burn	sSUB, Burn	Substrate and management
sSUB + daRad	sSUB, daRad	Substrate and radiation
sSUB + Precip	sSUB, Precip	Substrate and precipitation
sSUB + sArea	sSUB, sArea	Substrate and area
sSUB + sArea + Burn	sSUB, sArea, Burn	Substrate, area, and management
post hoc Additive	Burn, sSUB, m3Min, rTime, mEffort, daRad	post hoc
Manager Additive	Burn, sSUB, m3Min, rTime, mEffort	Manager variables

). 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² 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.

2.5. 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]:

$$p^{*}=1-{\left(1-p\right)}^k,$$

(1)

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:

$$k=\frac{\ln \left(-p^{*}+1\right)}{\ln \left(1-p\right)}$$

(2)

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 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.

3. Results

3.1. 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.

3.2. 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. Akaike information criterion (AIC) results for the 27 binary mixed-effect logistic regressions for detection probability models. Data were collected from 713 search events on the Eastern Massasauga at Carlyle Lake, Clinton County, Illinois from 1999–2014.

Rank	Model	K	−2LL	AICC	ΔAICC	ωi	r2marg	r2cond
1	post hoc	13	−399.65	825.82	0.00	0.83	0.28	0.40
2	Manager	11	−403.31	828.99	3.17	0.17	0.27	0.39
3	Global	30	−389.07	840.87	15.05	0.00	0.33	0.47
4	tEffort + Burn	5	−416.62	843.33	17.51	0.00	0.25	0.34
5	tEffort + sArea + Burn	6	−415.90	843.92	18.10	0.00	0.25	0.35
6	SearchEffortAdd	6	−418.55	849.21	23.39	0.00	0.23	0.35
7	tEffort + sArea	5	−420.73	851.55	25.72	0.00	0.22	0.33
8	Calander	7	−444.23	902.62	76.79	0.00	0.11	0.33
9	sSUB + sArea + Burn	7	−450.72	915.61	89.78	0.00	0.08	0.31
10	sSUB + Burn	6	−452.62	917.36	91.54	0.00	0.06	0.29
11	m3Min + daRad	6	−453.60	919.32	93.50	0.00	0.03	0.28
12	sSUB + daRad	7	−452.58	919.33	93.50	0.00	0.04	0.31
13	sSUB + Precip	6	−453.84	919.81	93.98	0.00	0.03	0.29
14	3DayTemps	9	−450.85	919.96	94.14	0.00	0.04	0.29
15	m3Min + Precip	7	−452.94	920.03	94.21	0.00	0.04	0.29
16	sSUB + sArea	6	−454.19	920.50	94.68	0.00	0.03	0.29
17	m3Mean + daRad	7	−454.20	922.56	96.73	0.00	0.03	0.28
18	pMin + Precip	6	−455.34	922.79	96.97	0.00	0.02	0.26
19	ProximateTemps	7	−454.74	923.65	97.82	0.00	0.03	0.28
20	m3Mean + Precip	6	−455.87	923.86	98.03	0.00	0.02	0.26
21	pMin + daRad	7	−455.15	924.45	98.63	0.00	0.03	0.27
22	Burn	4	−458.28	924.62	98.80	0.00	0.03	0.25
23	pMean + Precip	6	−456.38	924.87	99.05	0.00	0.02	0.25
24	PDayTemps	9	−453.71	925.67	99.85	0.00	0.03	0.27
25	pMean + daRad	7	−456.24	926.64	100.82	0.00	0.02	0.26
26	Null	3	−461.16	928.36	102.53	0.00	0.00	0.24
27	EvapCooling	10	−457.74	935.79	109.96	0.00	0.01	0.26

). Variables included in the model were Burn, daRAD, mEffort, rTime, sSUB, m3Min, and all associated quadratic terms (

Table 4. Summary statistics from the five variables which comprised the top post hoc model using data from search events on the Eastern Massasauga at Carlyle Lake, Clinton County, Illinois from 2002–2014.

Variable	Mean	SD	Min	Max	95% Confidence Interval
					Lower	Upper
sSUB	18.12	4.52	−0.10	34.00	17.79	18.45
daRad	806.34	121.91	144.00	1014.00	797.35	815.25
m3Min	5.84	4.98	−7.20	19.80	5.46	6.19
rTime	0.56	0.08	0.30	0.79	0.56	0.57
mEffort	53.19	37.38	5.00	224.00	50.41	55.90

). The post hoc model explained a high amount of variance in detection probabilities (r²_cond = 0.40; Table 3) with most of the variance explained by the fixed effects (r²_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²_cond = 0.39; Table 3) mostly by the fixed effects (r²_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²_cond = 0.47; Table 3), again, with most of the variance explained by the fixed effects (r²_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², rTime ≈ midday, sSUB ≈ 18 °C, m3Min ≈ 6 °C, and mEffort ≈ 50 min (Table 4; Figure 2).

Whether a site received a burn or not was significant, but prediction intervals (PIs) broadly overlapped. Burned sites had a 0.60 detection probability (95% PI = 0.46, 0.72) and unburned sites had a 0.44 probability (95% PI = 0.30, 0.59; Figure 3). Thus, burning somewhat aided the detection of Eastern Massasaugas.

When considering sSUB, the PIs were wider at higher values (due to the lack of surveys over ~30 °C) but narrower at cooler temperatures despite having few surveys under those conditions (

Table 5. Standardized parameter estimates, standard errors (SE), and confidence intervals (CI) from the top post hoc and Manager binary mixed-effect logistic regressions for detection probability. Data were collected from 713 search events on the Eastern Massasauga at Carlyle Lake, Clinton County, Illinois from 2002–2014. The best parameter model includes only terms with confidence intervals not bounding zero from the global model.

Parameter	post hoc 95% CI				Manager 95% CI
	Estimate	SE	Lower	Upper	Estimate	SE	Lower	Upper
Intercept	0.40	0.28	−0.15	0.96	0.31	0.27	−0.22	0.85
Burn	−0.65	0.28	−1.19	−0.10	−0.65	0.28	−1.19	−0.11
sSUB	0.47	0.12	0.22	0.71	0.31	0.27	−0.22	0.85
sSUB2	−0.20	0.07	−0.34	−0.06	−0.65	0.28	−1.19	−0.11
daRad	−0.36	0.13	−0.62	−0.10	-	-	-	-
daRad2	−0.11	0.07	−0.25	0.03	-	-	-	-
m3Min	0.11	0.11	−0.11	0.32	0.31	0.27	−0.22	0.85
m3Min2	0.12	0.09	−0.05	0.30	−0.65	0.28	−1.19	−0.11
rTime	−0.14	0.11	−0.35	0.07	0.31	0.27	−0.22	0.85
rTime2	−0.16	0.08	−0.32	0.00	−0.65	0.28	−1.19	−0.11
mEffort	1.07	0.13	0.81	1.33	0.31	0.27	−0.22	0.85

; Figure 1 and Figure 2). The general pattern illustrated higher detection probabilities between 18 and 26 °C (Figure 2). At the lowest sSUB (6.8 °C), detection probabilities were 0.10 (95% PI = 0.03, 0.27) but attained a maximum of 0.59 (95% PI = 0.46, 0.72) at 23.4 °C sSUB (Figure 2).

For daRad, the PIs were wide at lower values (due to the lack of surveys at conditions under ~700 J/s·m²), but the general pattern was a decrease in detection probability with increasing daRad (Table 5; Figure 1 and Figure 2). Detection probabilities for daRad peaked at 0.60 (95% PI = 0.45, 0.73) under ~620 J/s·m² (Figure 2). Then, as daRad increased to 1013 J/s·m², detection probabilities decreased to 0.31 (95% PI = 0.16, 0.50).

As m3Min increased, so did detection probabilities (Figure 3). However, at the extremes of m3Mins, we had greater uncertainty due to fewer surveys above ~18 °C and below −1 °C (Figure 2). With m3Min at −5.2 °C, we observed relatively high detection probabilities of 0.61 but had little predictive confidence in them (95% P.I. = 0.37, 0.82; Figure 3). We observed our lowest detection probabilities of 0.52 (95% P.I. = 0.40, 0.64) at an m3Min of 3.6 °C (Figure 3). Detection probabilities then steadily increased and peaked at 0.80 (95% PI = 0.50, 0.94) when m3Min was 19.8 °C (Figure 3).

We observed a peak in detection probability around midday, but again had greater uncertainty at earlier and later rTimes due to sampling (Figure 2 and Figure 3). With our earliest surveys at 0716 h, we had a detection probability of 0.25 (95% PI = 0.06, 0.64; Figure 3). However, detection probabilities were greater than 0.50 between 1042 and 1433 h, with the peak detection of 0.53 (95% P.I. = 0.41, 0.66) occurring at 1237 h (Figure 3). Thereafter, detection probabilities decreased to 0.17 (95% PI = 0.05, 0.46) at 1856 h (Figure 3).

Considering the mean effort per searcher, we observed a sigmoidal increase where detection probabilities had high predictive confidence along the full range of estimates and nearly achieved 1.00 (Figure 2 and Figure 3). At the lowest end of mEffort, 5 min resulted in a 0.22 (95% PI = 0.14, 0.33) detection probability (Figure 3). As mEffort increased, detection probability rapidly increased to 0.50 (95% PI = 0.37, 0.62) with 49 min of mEffort, and a 0.95 (95% PI = 0.87, 0.98) detection probability was achieved with 149 min (Figure 3). Additional mEffort showed diminishing returns, whereby a 0.99 (95% PI = 0.95, 1.00) probability was achieved at a 196 min and maxed out at 0.99 (95% PI = 0.97, 1.00) at 224 min (Figure 3).

3.3. 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⁻⁵ to 2.54 × 10⁻²⁸ by year (

Table 6. Overall probabilities and chances that Eastern Massasaugas were not detected during surveys at a historic site in the Chicago region, Illinois, United States of America (USA) (Table S3, Supplementary Materials), assuming the site was occupied.

Year	Probability Missed but Present			Odds Missed but Present
	Mean	Lower	Upper	Mean	Lower	Upper
2008	2.54 × 10−28	1.29 × 10−16	1.46 × 10−50	1:3.94 × 1027	1:7.77 × 1015	1:6.84 × 1049
2009	2.56 × 10−5	4.40 × 10−4	7.92 × 10−11	1:3.91 × 104	1:2.27 × 103	1:1.26 × 1010
2010	5.23 × 10−6	1.74 × 10−3	6.88 × 10−9	1:1.91 × 105	1:5.73 × 102	1:1.45 × 108
2015	1.80 × 10−8	3.95 × 10−4	1.46 × 10−15	1:5.57 × 107	1:2.53 × 103	1:6.84 × 1014
Overall	6.09 × 10−46	3.90 × 10−26	1.17 × 10−83	1:1.64 × 1045	1:2.57 × 1025	1:8.58 × 1082

). When considering all surveys, the probability was 2.56 × 10⁻⁴⁶ (95% CI, 1.17 × 10⁻⁸³, 3.90 × 10⁻²⁶) or, roughly, an average of a one in 1.64 × 10⁴⁵ 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.

4. 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.

4.1. 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–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.

4.2. 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:30 and 16:00 h (0.43–0.50), between substrate temperatures of 16.7 and 30.2 °C (0.49–0.59), and between mean minimum three-day temperatures of 15.9 and 19.7 °C (0.70–0.80), under overall solar radiation of 373–844 J/s·m² (0.50–0.60), for 1.8–2.4 h (0.85–0.95) per searcher. Under the Manager model (excluding daRad), the combinations of variables we present above should result in mean detection probabilities from 0.88 to 0.96, but with greater uncertainty.

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.

5. 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.