Occupancy and Detection of Eastern Massasauga Rattlesnakes (Sistrurus catenatus): Implications for Evaluating Population Recovery Efforts

There is a need to improve conservation translocation techniques for reptiles (Choquette et al. 2023), including evaluating appropriate monitoring methods to gauge translocation success (Germano et al. 2014). Effective posttranslocation monitoring, although challenging (Berger-Tal et al. 2020), can provide valuable information on population trajectories (Nicoll et al. 2021). Standardized estimates of a species’ distribution or abundance before and after translocations are essential, because these can inform whether a project was effective at long-term population establishment and persistence (IUCN/SSC 2013), assuming evaluations are made a sufficient amount of time after translocations (Seddon et al. 2012). Although reintroduced populations can grow initially through increasing density or increasing range (Armstrong and Seddon 2008), monitoring changes in abundance of certain taxa is inhibited by low detection rates (Steen 2010; Durso et al. 2011). For populations of cryptic species, monitoring changes in distribution (e.g., Bailey et al. 2004; Sewell et al. 2012) may be a more appropriate method for evaluating long-term translocation outcomes. Occupancy modeling is one approach that could provide distribution estimates both before and after translocations, to demonstrate whether a statistically significant increase in occupancy has occurred.

Occupancy modeling has grown in popularity as a robust and effective means to estimate the distribution (i.e., occupancy) of cryptic animals at the landscape scale, including snakes, because it accounts for imperfect detection (Durso et al. 2011; MacKenzie et al. 2017). Repeated standardized surveys at replicate sites allow practitioners to assign a probability of occupancy for each site based on the number of surveys conducted and the average detection probability (P) of the target species. To ensure unbiased occupancy estimates, five main model assumptions must be met: sites are closed to changes in occupancy, detections are independent, species are correctly identified, and there is no unmodeled heterogeneity in either occupancy probability (psi) or P (MacKenzie et al. 2017). Furthermore, the probability of detection should be ≥0.15 to avoid biased occupancy estimates generated when sites with low P are indistinguishable from sites where the species is truly absent (Bailey et al. 2004). Probability of detection in squamate reptiles may be influenced by many factors, including population density (Kéry 2002; Durso et al. 2011), survey year (Seddon et al. 2011), temperature (Harvey 2005), survey effort or method (Bartman et al. 2016; Bradke et al. 2018; Lardner et al. 2019), and time of year (Kéry 2002; Sewell et al. 2012). These factors complicate the development of an appropriate monitoring regime, particularly when a species is difficult to find.

Eastern Massasaugas (Sistrurus catenatus) are globally imperiled (G3; NatureServe 2022), small and cryptic pit-vipers that can be difficult to detect relative to sympatric snake species (Perelman et al. 2022). In an effort to increase survey efficiency, past protocols provided recommendations on the most suitable survey conditions, time frames, methods, and vegetation communities for detection and provided minimum survey effort targets (Casper et al. 2001; OMNRF 2016a). These protocols did not, however, stipulate an ideal number of searchers during surveys, describe seasonal variations in detection rate, or agree on whether visual encounter surveys ought to be supplemented with artificial cover object surveys (e.g., Bartman et al. 2016).

Recent studies of Massasauga P and occupancy in Illinois and Michigan complemented prior survey protocols by identifying variables affecting Massasauga P (e.g., Julian date, search effort, and substrate and air temperature; Shaffer et al. 2019; Crawford et al. 2020; Thacker 2020; Thacker et al. 2023); however, knowledge gaps remain. For example, Massasauga detection probabilities have not been estimated using data from the entire active season (spring, summer, and fall), or after accounting for a trap response. Snakes that are repeatedly found by surveyors in the same locations (i.e., “trap happy”; Bradke et al. 2018) may result in violating the independence assumption of occupancy modeling (MacKenzie et al. 2017) and inflating P estimates. Also, Massasauga detection probabilities estimated in “medium” or “large” subpopulations (e.g., Faust et al. 2011; Crawford et al. 2020; Thacker 2020), or in quasi-controlled settings with telemetered snakes (Shaffer et al. 2019), may not be transferable to small subpopulations at low densities (Kéry 2002; Royle and Nichols 2003). Experts have estimated that more than half of known Massasauga subpopulations are “small” (≤25–50 females; Faust et al. 2011), and a small subpopulation is presumably more consistent with a posttranslocation scenario. A comparison of results from standardized surveys conducted across all seasons, with a variable number of surveyors, and using different survey methods, while accounting for the effect of a trap response, would provide quantitative evidence to support the refinement of current rattlesnake monitoring protocols. Valuable conservation dollars could thus be directed toward the most efficient search techniques, potentially reducing material and human resource costs, increasing confidence in survey results for sites without detection, and shortening the time required to evaluate snake translocations.

Our goal was to estimate P and occupancy with a high level of accuracy and precision, in a small and endangered subpopulation of Eastern Massasaugas before experimenting with conservation translocations as a recovery tool, to determine the most appropriate monitoring regime. Because of a high risk of local extinction and low viability (COSEWIC 2012), recovery strategies urged further research into the effectiveness of conservation translocations as a management tool for increasing abundance, distribution and viability of small subpopulations of Massasaugas (Parks Canada Agency 2015; OMNRF 2016b). Our research questions were (1) which survey season (spring, summer, or fall), survey method (visual encounter or artificial cover object), and number of surveyors per survey (single or team) offer the highest P, while accounting for the effect of trap response; and (2) would the current sampling regime be adequate for monitoring changes in occupancy of this species after translocation?

Over 6 yr, we conducted repeated standardized surveys at 40 sites by using two methods, in three seasons, and with different numbers of surveyors, and we analyzed detection data using occupancy modeling. Surveys targeted Massasaugas within the Ojibway Prairie Complex and Greater Park Ecosystem, in Ontario, Canada (42.26°N, 83.07°W; datum = WGS84; Fig. 1), which represents the southernmost Canadian subpopulation. The federally and provincially endangered Ojibway Prairie Massasauga subpopulation is ecologically, genetically, and geographically distinct in Canada, yet before our study, the range of this small subpopulation had declined by 82% in three decades, and abundance was critically low (≤40 mature snakes; Faust et al. 2011; COSEWIC 2012; COSSARO 2013). Conservation translocations were conducted in the 1990s and 2000s in an effort to augment the local population, but were presumed unsuccessful (19 snakes released between 1990 and 1994 [Pither 2003]; 27 snakes released in 2006 [Harvey et al. 2014]). More recently, Lincoln–Petersen estimates generated an average abundance of only 10 (95% confidence interval; [CI] = 6–80) remaining adults and subadults across ca. 20 ha in any given 2-yr period during a long-term study (J.D. Choquette and C. Fournier, personal observations).

Survey sites (n = 40) were identified as areas or discrete patches of open-canopy or semiopen-canopy vegetation that were 2.0–2.5 ha (average, 2.4 ha) and found within the recent historical subpopulation range of Massasaugas at the study area (Fig. 1). A habitat suitability model (Choquette et al. 2020), a Massasauga Critical Habitat map (Parks Canada Agency 2015), aerial photos (LIO 2006), and ground-truthing were used in combination to guide the site selection process. Most survey sites (83%; n = 33) were within or directly adjacent to Critical Habitat and we assumed these sites were representative of all Critical Habitat because we selected the majority of available sites therein and followed a probabilistic sampling scheme (MacKenzie and Royle 2005; Supplemental Methods S1, available online). Open and semiopen vegetation communities were targeted (e.g., meadows, prairies, thickets, and savannahs) because Massasaugas preferentially use vegetation types with relatively low canopy cover during the active season (COSEWIC 2012). Size of selected sites approximated the average size of the remaining patches of open or semiopen vegetation in the study area and was similar to the mean annual Massasauga home range size in nearby southern Michigan, USA (e.g., 2.5 ha; Bissel 2006). A geographic information system was used to map boundaries and estimate area of each survey site (ArcGIS 9.2, ESRI, Redlands, CA, USA). Site boundaries were generally discrete (e.g., edges of roads, trails, or forests), except where the presence of relatively extensive patches of open-canopy vegetation at one of the Critical Habitat blocks necessitated the delineation of arbitrary site boundaries (Supplemental Methods S1). All sites were similar in size (within 0.5 ha) to reduce impacts of survey area on occupancy or P (Thacker 2020).

Surveys occurred over a 6-yr period, and a subset of all sites was repeatedly surveyed in any given year. Site 9, for example, was surveyed 30 times overall, including 9 times in 2014, 14 times in 2015, and 7 times in 2016. We aimed to complete 19–29 surveys at each site over the 6 yr, of which 4–6 surveys were to be completed in each of two summers to account for the biennial reproductive cycle of Massasaugas (COSEWIC 2012; Supplemental Methods S2, available online). Relatively greater survey effort was invested at three specific sites across all 6 yr to produce repeat Massasauga captures for demographic estimates and to determine factors associated with increased P. These three sites were chosen based on their proximity to the most recently confirmed observations, and surveys were conducted therein during the spring egress, summer gestation, and fall ingress periods (Supplemental Methods S3, available online). At all other sites, surveys occurred mostly during the summer gestation period because (1) gravid females would have dispersed to gestation sites and have a relatively higher detection rate than other age classes (Hileman et al. 2018) and (2) any snakes that had been hibernating in closed-canopy vegetation outside of survey sites would have had time to migrate into patches of open-canopy vegetation. It was assumed that no colonizations or extirpations occurred at the site level during the 6 yr due to the strong site fidelity demonstrated in this species (COSEWIC 2012) and because the occupancy status of the three sites with greater survey effort remained constant during that period (note that an incorrect closure assumption may result in an overestimation of psi; MacKenzie et al. 2017).

Survey start and end dates in each calendar year varied due to staffing logistics and arrival of favorable weather conditions: they generally began in April and ended in September (median start and end date, 21 April and 22 September, respectively; start range, 9 April–2 May; end range, 23 August–5 October). All surveys took place within the documented activity period for the Massasauga in southeastern Michigan and southwestern Ontario (Pratt et al. 1993; Bissell 2006; Moore and Gillingham 2006; Bailey 2010). Thus, we presumed snakes were available for detection at occupied sites from the start to end of each survey period in any given year. In all years, regardless of season, each survey was generally conducted by a single surveyor on an infrequent basis (i.e., one to three surveys per site per week). From 2015 to 2018, however, spring surveys were generally conducted by a team of two or more surveyors at the same three sites more frequently (i.e., three to five surveys per site per week). Survey team composition varied annually, and at least one individual trained in our protocols was present on 100% of surveys across the study period.

Survey methods included time-constrained standardized visual encounter surveys and artificial cover object surveys. Visual encounter surveys were 1.5–2.0 h in duration and were conducted following recommended survey protocols (Supplemental Methods S4, available online). Artificial cover object surveys were conducted during visual encounter surveys in all years except 2013 (Supplemental Methods S5, available online). Survey sites were covered as evenly and thoroughly as possible, although portions with suitable microhabitats such as basking features and animal burrows were given relatively greater attention. Survey start times (AM vs. PM), surveyor(s), and start/end points were rotated to reduce heterogeneity in P arising from time of day, surveyor, and direction of travel. We attempted to meet occupancy modeling assumptions by standardizing our study design and addressed the remainder of assumptions using model covariates (Supplemental Methods S6, available online).

Detailed observation and weather data were recorded in the field generally following recommendations in a 2001 survey protocol (Casper et al. 2001). Surveys were conducted during a predefined range of recommended weather conditions (Supplemental Methods S4) to reduce the impact of associated variables on detection as much as possible. Shaded air temperate was recorded at waist height to the nearest degree Celsius with a digital thermometer (model 317C, Test Products International Inc.; accuracy, ±1°C; model 9835, Taylor Precision Products Inc.; accuracy, ±0.6°C), wind speed was coded using the Beaufort scale, and cloud cover was recorded via subjective assignment of a weather category (e.g., sunny, overcast, etc.) and by estimating cloud cover to the nearest 10%. All weather data were recorded immediately before and after each survey. Any missing air temperature or weather condition data were retrieved using Environment Canada’s historical hourly weather data for Windsor International Airport (∼7.5 km from study site) for the hour nearest the time of snake observation or survey start/end.

Occupancy modeling analyses were completed in program Presence (v2.13.18, released December 2021; USGS Patuxent Wildlife Research Center, Blacksburg, VA, USA). Only results from standardized visual encounter surveys were included in analyses. Results from 237 nonstandardized surveys and site visits were excluded from analyses, even if the target species was detected, because they did not adhere to our standardized survey protocol. For example, we excluded surveys that were conducted with air temperatures or wind speeds outside of the protocol based on the average of start and end measurements, repeat surveys at the same site on the same day, site visits as part of related conservation projects, and surveys that were <1.5 h long. Any surveys that were >2.0 h long were included in analyses based on results from the first 2.0 h only. Massasauga observations that did not occur during standardized visual encounter surveys (including those reported to us by park users) were labeled “incidental.” Six surveys ended early after a Massasauga was detected (i.e., before 1.5 survey hours), but these were so few that their inclusion had a negligible impact on P estimates (n = 6; 0.3% of standardized surveys; increase in P of 0.006).

The detection probabilities were first estimated separately for each calendar year of study by using a single season analysis, 1000 bootstraps, no covariates on occupancy or detection (i.e., a null model), and a detection history (DH) wherein sampling occasions were based on calendar day. Estimates for single years were generated to determine the optimum number of repeat surveys (Supplemental Methods S2) and served as null models against which we compared covariate-fitted models. Regardless of year, 1–11 (mean, 3) sites were surveyed during any given sampling occasion (i.e., calendar day) and results for sites not surveyed during a sampling occasion were considered “missing” observations. Yearly variations occurred in the total number of sites surveyed (mean, 19; range, 9–31), sampling occasions (mean, 102; range, 79–114), and missing observations (mean, 1616; range, 785–2675). The daily DHs from each year were combined to create a complete 6-yr dataset that was then reformatted such that each sampling occasion represented the first to the nth survey at each site following the chronological sequence of surveys (Supplemental Methods S7, available online). In other words, the first survey at each site was considered sampling occasion 1, the second survey was sampling occasion 2, and so on, regardless of the calendar day when the surveys were conducted.

For the combined 6-yr dataset, we conducted separate single season analyses using thee different DHs (daily, weekly and monthly; Table 1), but all were based on the same underlying detection data. The daily DH is as described above, whereas the weekly and monthly DHs were coded based on the pooled weekly or monthly results of all surveys at each site. This allowed for a comparison of results at different temporal scales, each providing a unique tradeoff between the number of missing observations and the temporal precision of the data (Table 1). One site covariate and six survey covariates were selected a priori based on their suspected influence on psi or P, and all or a subset of these were included in relevant model sets (Supplemental Methods S3). For the site covariate, general location, each site was assigned to one of four categories based on its proximity to a distinct Critical Habitat block (Critical Habitat Block 1, 2, 3, and other). The survey covariates were generally categorical and as follows: season (meteorological spring, summer, or fall), survey frequency (daily DH only, <5 or ≥5 d since prior survey of a site), calendar year, cloud cover (daily DH only, sunny or cloudy), and search effort (daily DH, single or team surveys; weekly DH, one or two or more surveys per week; monthly DH, continuous based on number of surveys per month). We also included a covariate on trap response, based on whether a trap-happy snake was present or absent at a site during a survey. Trap-happy snakes were individuals found within 10 m of the same location on at least three successive occasions (Supplemental Methods S3). Additive models using logical combinations of these base covariates were also used in model sets. To avoid model overparameterization, and because we followed a standardized survey protocol regarding weather conditions (Supplemental Methods S4), and surveyed sites of similar size and with suitable open-canopy vegetation, we excluded time of day, air temperature, wind speed, site size, and vegetation community type as model parameters. A unique global or subglobal model was used in each analysis to assess model fit using 1000 bootstraps (Lamothe et al. 2021), and from this, c-hat was adjusted to account for model overdispersion. The Quasi-Akaike’s Information Criterion (QAIC) was used to evaluate models.

A preliminary step was taken before running the daily, weekly, and monthly analyses to identify relative importance of the six survey covariates and eliminate those with low support. This reduced the eventual number of parameters being estimated, thereby avoiding model overparameterization (MacKenzie 2006; MacKenzie et al. 2017). Using a general model for occupancy psi(location), a suite of models was run using all relevant combinations of single and paired covariates on detection (i.e., for the daily DH each of the six covariates on detection were represented in six models). Model weights were summed for the group of models represented by each covariate, and covariates with relatively low support (e.g., summed model weight ≤0.001) were excluded from the next step in the analyses. The final daily, weekly and monthly analyses included one site covariate and combinations of three or four survey covariates, as well as the null models. The most general (global) model was included in all analyses. Overparameterized models were removed from each model set, which occurred when the number of significant digits was less than six (the default in Presence), or when negative SEs appeared in the variance–covariance matrix. The result was a set of 11, 9, and 22 models in the daily, weekly, and monthly analyses, respectively. Results were presented for all models with substantial support (ΔAIC = 0–2), the next highest model (ΔAIC > 2, if applicable), and the null model. A critical appraisal was conducted to further explore the influence of survey covariates on P (Supplemental Results S8, available online). Models with no support (AIC > 10) were not presented (Supplemental Results S9, available online).

In total, 1742 standardized visual encounter surveys were completed across 40 sites, representing 3822 person-hours (Table 2). The median number of surveys per site overall was 29 (range, 5–255); however, this value was higher for the seven sites at Critical Habitat Block 1 (118; range, 30–255) and lower for the eight sites at Critical Habitat Block 2 (10; range, 5–29). Most sites (80%; n = 32) received the target of 19–29 standardized surveys overall, and the majority of sites (90%; n = 36) received the target of four to six summer surveys in at least two different calendar years (Supplemental Methods S2).

At least one Massasauga was detected during 150 out of the 1,742 standardized visual encounter surveys (8.6% of surveys), and all detections during visual encounter surveys occurred entirely within four sites, all within Critical Habitat Block 1 (naïve occupancy = 0.10; 4/40 sites). Detection rate at the four occupied sites was 18% (150/842 surveys). Overall, surveys were conducted most often by a single surveyor (74%), on an infrequent basis (63%), or in the summer (61%; Supplemental Results S10, available online). At the four occupied sites, a comparable proportion of surveys were conducted in spring (48%) and summer (45%), at the expense of fall surveys (7%), and single surveyors completed a similar proportion of surveys (49%) as team surveyors (51%; Supplemental Results S10). A trap-happy snake was only present during 15% of surveys at occupied sites (Table 2; Supplemental Results S10). No Massasaugas were detected within Critical Habitat Block 3, where conservation translocations were previously conducted (Pither 2003; Harvey et al. 2014).

We recorded 225 Massasauga observations, representing at least 39 unique individuals; 185 of these observations occurred during standardized visual encounter surveys, and >1 snake was found during some surveys. Forty of the observations were incidental (i.e., they did not occur during standardized surveys), but all occurred within, or very close to, the four occupied sites. No incidental observations were included in occupancy modeling analyses. In each year, and excluding neonates, more than half of observations were of individuals already marked in that year (range, 54–89%). Included in the grand total are two observations of dead adult Massasaugas (2013, 2018), two adult skin sheds (both in 2018), and five litters of neonates (2013, 2015, 2016, and 2017; representing 34 observations). Zero to two occupied sites were used annually by Massasaugas for reproduction, which was determined based on the detection of neonates near a postpartum female, and all four occupied sites were used for reproduction during the study.

Survey conditions were summarized below using data from 196 Massasauga observation records. For clarity, all 225 Massasauga observations were reduced to 196 after combining observations of multiple neonates found together into a single observation record for each litter. Observations occurred from 10 April to 12 September, with a median date of 1 June (i.e., 151st calendar day) and litters of neonates were observed from 16 to 30 August. A large proportion of observations occurred in May (40%), June (19%) and August (17%), compared with April (9%) and September (2%). Observations were made equally in the spring (49%) and summer (49%), whereas very few occurred in the fall (2%). Excluding neonates, a similar number of unique individuals were observed in the spring of each year (mean, 5) compared with the summer (mean, 4). Mean shaded air temperature for all Massasauga observations was 21°C (SD, 4.2; range, 10–33°C; Supplemental Methods S4), and the peak number of observations occurred between 22 and 25°C (n = 73; 38%). Median time of day for all observations was 1335 h (range, 0903–1940 h), with the peak between 1100 and 1159 h (n = 39, 20%). Mean cloud cover for all observations was 54% (SD, 33%; range, 0–100%), with the peak when cloud cover was 100% (n = 38; 20%). Time of day and/or date was not recorded for two incidental Massasauga observations.

The vast majority of Massasaugas were detected by sight (91%), with very few detected by sound (8%) or found under existing wood and metal rubbish (2%). During surveys, Massasaugas were only detected using visual encounter search methods, and none were found under our standardized artificial cover objects. We conducted 1343 cover object surveys, 615 of which occurred at sites occupied by Massasaugas, and these generated detections of five other snake species: Pantherophis vulpinus, Storeria dekayi, Storeria occipitomaculata, Thamnophis butleri, and Thamnophis sirtalis.

Annual estimates of P using a null model ranged from 0.05 to 0.35 and were above 0.15 in 2013, 2017, and 2018 (Table 2). The combined 2013–2018 dataset included 34–255 survey occasions and 974–8458 missing observations (Table 1). Sample covariates on search effort and trap response were of relatively high importance regardless of the version of DH used, whereas season was of high importance only with the monthly DH (Table 3). All other covariates on detection, which included cloud cover, survey frequency, and calendar year, were of relatively low importance (summed model weight ≤ 0.001; Table 3).

Results from the top model using the daily DH indicated that P was influenced by season with additive effects from search effort and trap response (Table 4). Estimates of P were significantly higher for surveys influenced by trap response vs. those without, for team surveys vs. single surveys, and for summer or spring surveys vs. fall surveys (i.e., 95% CIs did not overlap; Table 4; Supplemental Results S9). Furthermore, trap response appeared to have a stronger effect on P than season or search effort (Supplemental Results S8). The highest overall P occurred during the summer, with team surveys when a trap-happy snake was present (0.73, 95% CI = 0.60–0.83). When surveys were not influenced by trap response, team surveys achieved a P at or above the 0.15 threshold in spring and summer (0.15 and 0.26, respectively), but not in fall (0.08; Fig. 2).

Regardless of whether the weekly or monthly DH was used, trap response and search effort remained important covariates on P (Table 4). An additive effect between search effort and trap response using the weekly DH substantially improved model fit with only one additional parameter, as indicated by a decrease in the −2 log-likelihood (−2l) value of ca. 30 or 54, compared with weekly models with only one or the other covariate (Supplemental Results S9). Based on the top weekly model, when surveys were not influenced by trap response, P was ≥ 0.15 only when two or more surveys were completed per week (Table 4). The highest P overall, using the weekly DH, occurred during weeks when two or more surveys were completed and when a trap-happy snake was present during a survey (0.76, 95% CI = 0.63–0.86). Although the inclusion of the calendar year covariate did improve model fit based on the −2l value for both weekly and monthly DHs, it required five additional parameters and did not generate significantly different P estimates between years (Table 4). Using the second-best monthly model, and when surveys were not influenced by trap response, a significantly higher P occurred in the spring (0.59) and summer (0.39) than in the fall (0.04).

Based on the top model using the daily DH, the estimated Massasauga psi in the study area was 0.13 (95% CI = 0.05–0.30), which was very similar to psi estimates from the highest ranking models using the weekly DH (Table 4). Because we sampled a finite population of sites (MacKenzie et al. 2017), a psi of 0.13 estimated the proportion of our study sites that were occupied by Massasaugas (i.e., 5/40 sites). Note that the SEs estimated by program Presence, and associated confidence intervals (Table 4), are imprecise for finite populations (MacKenzie et al. 2017). Massasaugas were only detected at four sites (naïve occupancy, 0.10). Therefore, results suggest this species was present and remained undetected at a fifth site during the study period.

The site-specific probability of occupancy at sites where we did not detect Massasaugas was very low overall (mean psi-conditional, 0.03; range, 0.01–0.09; n = 36). When grouped by general location, and for sites with nondetection, site-specific probability of occupancy was highest at sites in Critical Habitat Block 2 (mean psi-conditional, 0.06; range, 0.02–0.16; n = 8). Furthermore, the four sites with the greatest probability of Massasauga occupancy (psi-conditional approximately, 0.09; 95% CI approximately, 0.03–0.22) were all situated in Critical Habitat Block 2 and were the sites surveyed the fewest number of times (five or six surveys each). Massasaugas were thus most likely to have escaped detection at a site in Critical Habitat Block 2. Assuming all five occupied sites were in Critical Habitat, and extrapolating to all available sites therein (Supplemental Methods S3), Massasaugas occupied an estimated 12% of Critical Habitat during our study.

An equally high-ranking model using the daily DH suggested that psi was not constant across the study landscape but was influenced by general location (the psi [loc] model; Table 4). This model differed from the top model by only one parameter and had a very similar −2l value (Δ−2l < 4); therefore, we deemed the general location parameter as uninformative when using the daily DH, and the second-ranking model as not competitive (Arnold 2010). The two top monthly models also included the general location covariate. Both models provided an estimated psi for sites in Critical Habitat Block 1 of 0.59 (95% CI = 0.23–0.86), whereas a psi of 0.00 (95% CI = 0.00–1.00) was generated for all other sites (Table 4). We deemed occupancy outputs from those models as uninformative given psi could not be estimated with any meaningful level of certainty for sites outside of Critical Habitat Block 1. Presumably, limitations were caused by the use of a categorical covariate to account for spatial correlation, as opposed to a continuous covariate (e.g., Euclidian distance; MacKenzie et al. 2017).

Our results provide support for a strong influence of trap-response on Massasauga detection probability. The presence of trap-happy snakes likely explains the relatively high P estimated separately for 2013, 2017, and 2018 when using the null model, because yearly differences in P were not apparent after accounting for trap response and using the combined 2013–2018 dataset. Although the highest estimates of P in our study were for surveys influenced by trap response, prior knowledge of the location of individual snakes violates occupancy modeling assumptions, so a comparison of survey techniques is only meaningful in the absence of those effects. After accounting for trap response, standardized visual encounter surveys conducted during spring and summer and by a team resulted in the highest P (0.15–0.26). In other words, surveyors detected the species during 2 or 3 out of every 10 surveys, on average, at sites occupied by Massasaugas and when trap-happy snakes were not present.

Results from the daily DH provided evidence for a significantly greater P during team surveys compared with single surveys, for all seasons combined and also during spring and summer seasons only (i.e., 95% CIs do not overlap). P also significantly increased as the number of surveys completed per week increased from one to two or more. These results are consistent with those of other studies wherein a positive relationship between Massasauga P and search effort was reported (Shaffer et al. 2019; Crawford et al. 2020; Thacker et al. 2023). Results from both the daily and monthly DH also provide evidence of a seasonal effect on P. Regardless of the DH used, surveys during the fall were significantly less productive than those in other seasons because average estimates of P were ≤0.08 in all cases and this result contradicts prior recommendations to survey during the entire Massasauga active season. Average P was higher during summer than spring using the daily DH, but this result was reversed when using the monthly DH, and in both cases, differences were nonsignificant. Further study may reveal calendar month to be a more informative covariate than spring or summer. For example, we recorded the greatest number of snake observations per standardized survey in one spring and one summer month (0.26 snakes/survey in each of May and June), which was higher than encounter rates in April (0.14), July (0.15), or August (0.19).

Our method for conducting artificial cover object surveys was ineffective at sampling Massasaugas. Perelman et al. (2022) also found no Massasaugas under artificial cover objects in Pennsylvania, USA, and Bartman et al. (2016) found gravid females using artificial cover objects in Michigan as basking surfaces only. Conversely, Massasaugas have been detected by other researchers under larger plywood and metal cover objects than those that we used (Amber et al. 2021; Yagi et al. 2023; T. Preney, personal observation). Further study on the most beneficial type and size of artificial cover object for sampling Massasaugas is warranted in an effort to improve P estimates.

We did not incorporate all variables recorded during surveys into our models as covariates to avoid overparameterization and because we standardized our site selection approach and our survey protocol. Regardless, additional insights into most suitable survey conditions can be gleaned from the summary of our Massasauga observation data, which may be used to refine survey protocols and improve P during future surveys: (1) the majority of observations (95% of 194 observations) were made from 0900 to 1800 h, with a peak from 1100 to 1300 h (36%). Similarly, Crawford et al. (2020) found that Massasauga P reached a maximum at ca. 1300 h; (2) the majority of observations were made during air temperatures ranging from 14 to 29°C (95% of 194 observations), with a peak occurring during temperatures ranging from 22 to 25°C (38%). Shaffer et al. (2019), however, found Massasauga P was greatest during air temperatures ranging from ca. 13 to 21°C; and (3) Massasaugas were observed in similar numbers under the full range of cloud cover scenarios (0–100% cloud cover), with a slight peak in observations during surveys with 100% cloud cover (20% of 193 observations). Similarly, in a study by Thacker et al. (2023), cloud cover was not a strong predictor of Massasauga P during spring and summer surveys.

Occupied sites were spatially clustered within Critical Habitat Block 1 only, yet our top daily model suggested no significant effect of the general location covariate on occupancy. The Massasauga psi estimated by that model (0.13) is less than half that estimated for Massasaugas across southern Michigan (Thacker et al. 2023) and similar to occupancy probabilities estimated for rare aquatic snakes in South Carolina, USA (Durso et al. 2011). Occupancy estimates generated by our top models using both the weekly and daily DHs were higher than those estimated by the null models; therefore, estimated occupancy would have been even lower had we not included covariates on detection. In effect, the estimated number of our study sites occupied by Massasaugas increased from four to five after accounting for detection heterogeneity.

Occupancy estimates had poor precision for the psi(.) model regardless of the DH used (SE [psi]/psi = 48%; see Bailey et al. 2004), and for our top daily model, the estimated CI was large relative to the proportion of sites occupied (upper 95% CI ca. 12 sites). Thus, a substantial increase in Massasauga distribution would need to be documented in a posttranslocation scenario before we could confidently demonstrate that a statistically significant increase in psi has occurred. Because we sampled a large proportion of the sites of interest, a finite population correction such as the delta method might improve the precision of the psi SEs. For example, the delta method may reduce the SE on psi estimates by up to 25–39% (MacKenzie et al. 2017). Additional investigations are warranted to assess the extent to which the precision of our psi estimate could be improved, with the aim of reducing the statistical threshold needed to detect a posttranslocation change in occupancy. Certainly, our results provide a liberal estimate of the proportion of sites occupied by Massasaugas during the study period.

The suitability of occupancy modeling for estimating distribution and changes therein in cryptic animals can be challenged by low detectability, resulting in poor or uninformative occupancy estimates (MacKenzie et al. 2002; Guillera-Arroita et al. 2014; Durso and Seigel 2015). Therefore, at occupied sites, the survey methods should be able to detect the target species at least 15% of the time (Bailey et al. 2004). Based on our results, the 15% threshold was only met during team visual encounter surveys in spring and summer by using the daily DH, or when two or more surveys were conducted per week by using the weekly DH. The P may be further improved by focusing survey effort within certain months (e.g., May, June, or August); however, to avoid biased occupancy estimates, surveys should be spread across both spring and summer. For example, in 2014, Massasaugas were not observed at two sites until summer (June and July, respectively), and in 2017, Massasaugas were detected at four sites in spring (May), yet only one of those sites yielded summer detections.

Massasauga occupancy surveys should be designed to detect evidence of reproduction if results are to be used to gauge whether or not population establishment has occurred (Choquette et al. 2023). Evidence of Massasauga reproduction could include (1) spring detections of yearlings that were born the previous summer, (2) detections of “suspected” gravid females in spring or summer, and (3) detections of neonates with a postpartum female in mid to late August (Bissell 2006). The latter would provide the strongest evidence for current year reproduction at a survey site; however, it may require intensive search effort at suspected or known gestation sites to confirm. For example, multiple surveys per week at gestation sites in August may be needed to successfully intercept postpartum females with their litters before dispersal. In this case, a removal sampling design applied in summer could help to focus resources toward a subset of sites (MacKenzie et al. 2017).

Declaring absence of cryptic vipers from a site with a high level of confidence may require substantial levels of search effort, particularly when population density is low (Kéry 2002). With a P ranging from 0.25 to 0.30, OMNRF (2016a) recommended conducting 10 surveys over two active seasons before concluding absence with 95% confidence, but that greater survey effort would be necessary for low density populations. Our study was focused on a small, low-density population of Massasaugas, and given an average of 25 repeat surveys at the 36 sites where the species was not detected, we report a mere 3% chance, on average, of Massasauga occupancy at any of those sites. Alternatively, and assuming a P of 0.21 (average of spring and summer P from daily DH), 14 negative visual encounter surveys conducted by a team of at least two surveyors would be required to declare absence from a site with 95% confidence (Fig. 2).

We did not include an exhaustive list of environmental variables as covariates in our analyses because we chose to conduct our surveys during conditions already deemed suitable for Massasauga detection. Regardless, our top models may have omitted finer scale influences of certain variables on Massasauga P (e.g., air temperature and time of day; Shaffer et al. 2019; Crawford et al. 2020), and we obviously cannot rule out the effect of variables for which we did not explicitly account. For example, although we selected sites with open and semiopen canopies, we did not quantify canopy cover and thus did not directly account for its effect on occupancy (Thacker et al. 2023). Nor did we record substrate temperature or solar irradiance, two variables known to affect Massasauga P (Crawford et al. 2020; Thacker 2020). Finally, the way in which we coded the trap response covariate may have resulted in overly conservative (low) P estimates after accounting for trap response; we treated all surveys influenced by a trap-happy snake in the same way, ignoring whether a non–trap-happy snake was also found during a survey (two or more individuals were detected during 13 surveys, excluding litters of neonates with mother). Regardless, a monitoring protocol for an endangered species based on a conservative P estimate would reduce the risk of declaring false absences.

We did not cross validate our model results with a data set from another study area with similar vegetation communities and Massasauga abundance and therefore cannot fully comment on the transferability of our results (e.g., Wenger and Olden 2012). Some insights into potential transferability, however, were gained via critical analysis of model predictions using two subsets of our dataset a posteriori. Those results suggested that the effect of trap response was the most important variable influencing P, compared with the effect of season (spring vs. summer) or search effort (single vs. team surveys in the spring). Therefore, a trap response effect on P is presumed to be the most transferable of our results to other study systems.

Our study generated baseline distribution data that can be used to evaluate changes in the distribution of a small subpopulation of Massasaugas after population management activities such as conservation translocations. We found no evidence that Massasauga translocations previously conducted at one Critical Habitat Block 12–24 yr before our study were successful at population establishment; therefore, more intensive efforts are warranted. Results can also provide guidance on improving survey protocols to increase P for this and other cryptic species. Based on our site size and search time, spring and summer surveys conducted by at least two people, and at least twice weekly are recommended. Focusing surveys in May (spring) and June (summer) may increase detection rate, and surveys in August to detect neonates are recommended if aiming to confirm population establishment. Surveyors should avoid using the artificial cover object survey method of as presented herein and should avoid conducting visual encounter surveys during the latter part of the Massasauga active season (September and October), to avoid wasting effort when P is lowest. We stress the importance of explicitly accounting for the effect of trap response when planning occupancy studies and when modeling datasets; trap-happy snakes artificially inflate P, leading to a violation of occupancy modeling assumptions and the potential to misguide practitioners regarding the true amount of search effort required to detect small subpopulations of Massasaugas and potentially other viperids.

Supplemental material associated with this article can be found online at https://doi.org/10.1655/Herpetologica-D-23-00055.S1, https://doi.org/10.1655/Herpetologica-D-23-00055.S2, https://doi.org/10.1655/Herpetologica-D-23-00055.S3, https://doi.org/10.1655/Herpetologica-D-23-00055.S4, https://doi.org/10.1655/Herpetologica-D-23-00055.S5, https://doi.org/10.1655/Herpetologica-D-23-00055.S6, https://doi.org/10.1655/Herpetologica-D-23-00055.S7, https://doi.org/10.1655/Herpetologica-D-23-00055.S8, https://doi.org/10.1655/Herpetologica-D-23-00055.S9, https://doi.org/10.1655/Herpetologica-D-23-00055.S10.

Funding for this study was provided by a Natural Sciences and Engineering Research Council of Canada CREATE grant to Laurentian University (ReNewZoo 481954-2016), a Carolinian Canada Coalition grant to Wildlife Preservation Canada (WPC; CSI Internship, 2013), Ontario Ministry of Natural Resources grants to WPC (SARSF 119-13-WPC4 and 115-14-WPC3), Ontario Ministry of Natural Resources and Forestry grants to WPC (SAR-00097 and 109_18_WPC3), a Colleges and Institutes Canada grant to WPC (Clean Tech Internship 2016), Employment and Social Development Canada grants to WPC (Canada Summer Jobs internship in 2017 and 2018), Environment and Climate Change Canada grants to WPC (HSP6824WPC, HSP6988WPC, and HSP8062WPC), and a Ganawenim Meshkiki and Eastern Georgian Bay Initiative grant to WPC (EGBI 2021-08). No funders played a direct role in the design, development or undertaking of this study. This project has received funding support from the Government of Ontario; such support does not indicate endorsement by the Government of Ontario of the contents of this publication. We thank the numerous WPC field technicians and interns who assisted in completing the surveys and collecting the data used in this study: M. Bagnall, J. Barden, J. Dertinger, D. Dunlop, B. Eisner, L. Gagnon, K. Hall, E. Jolin, T. Lourie, P. Nutley, E. Primeau, and L. Valliant. The Ojibway Nature Centre, City of Windsor (K. Cedar and T. Preney) provided data from 11 standardized visual encounter surveys completed by their staff in 2017. We are also thankful to the many others who volunteered their time to set up survey sites and to assist with surveys, especially K. Cavanaugh, J. Hartley, R. Jones, S. Marks, and M. Virtanen. C. Fournier assisted with manuscript formatting, figures, and French translation. Advice on study design or interpretation of results was provided at different stages of this project by J. Hines, K. Lamothe, D. MacKenzie, and D. Sewell. Advice and supervision were also provided by: R. Heide, A. Lentini, A. Shulte-Hostedde, J. Steiner, E. Williams, L. Woolaver, Jr., A. Yagi, and members of the Ojibway Prairie Reptile Recovery Working Group. Four anonymous reviewers provided feedback on an earlier version of this manuscript. City of Windsor, Province of Ontario, Town of LaSalle, and private landowners provided land access. Work with Massasaugas was conducted with required permits and authorizations under the Ontario Endangered Species Act (AY-B-010-13, M-102-1361488321-14, M-102-1797287827-15, M-102-5125809052-16, M-102-1175396436-17, M-102-2234438415-18), the Ontario Fish and Wildlife Conservation Act (1073890-13, 1076705-14, 1076705-15, 1082420-16, 1085806-17, 1089102-18), and the Ontario Ministry of Natural Resources & Forestry Wildlife Animal Care Committee (299-13, 299-14, 299-15, 299-16, 299-17, 299-18).