Snakes inhabit a range of climates, and those living in temperate ecosystems must adapt behaviorally and physiologically to the seasonal changes in their environment (Gregory, 1982; Huey, 1982; Ultsch, 1989). Temperate snakes have adapted to their climate by employing a biphasic life-history strategy (Huey, 1982). During the above-ground or the active season, surface temperatures are conducive for movement, foraging, growth, and reproduction (Harvey and Weatherhead, 2006a). When surface temperatures are suboptimal, snakes retreat below ground during the overwintering or hibernation phase (Cowles, 1941; Gregory, 1982; Ultsch, 1989; Costanzo and Lee, 2013). In some latitudes or altitudes, the overwintering phase may account for more than half of the animal's annual life cycle (Gregory, 1982). A typical snake overwintering site consists of a subterranean space with access to the surface that includes crayfish and mammal burrows, rotting tree roots, fractures in bedrock, caves, or holes weathered by wind, or water (Weatherhead and Prior, 1992; Johnson, 1995; Kingsbury and Coppola, 2000; Harvey and Weatherhead, 2006b). The depth, width, or ecological characteristics of the subsurface wintering site are not well understood (Prior, 1997), however, likely vary with soil type, geology, hydrology, vegetation, climate, and topography.

Radiotelemetry has provided insights into hibernation thermal ecology, the locations of hibernacula, site fidelity (Macartney et al., 1989; Johnson, 1995; Sage, 2005; Harvey and Weatherhead, 2006b), and the timing of ingress and egress (Nordberg and Cobb, 2016). Surgically implanted temperature data loggers or hygrochron data loggers fused to rattles and doubly labeled water have also been used to study the ecophysiology of hibernation (Agugliaro, 2011; Nordberg and Cobb, 2016). Because small numbers of relatively old and large snakes are monitored during radiotelemetry studies (Johnson, 1995; Harvey and Weatherhead, 2006b; Refsnider et al., 2012; Zappalorti et al., 2014), survival rates are difficult to determine (Jones et al., 2012). Radiotelemetry studies cannot locate all hibernation sites used by solitary hibernators. Therefore, estimating population level impacts using radiotelemetry alone may be inaccurate (Jones et al., 2012). Furthermore, the invasive surgical procedures associated with radiotelemetry itself may contribute to snake mortality (Lentini et al., 2011). Locating hibernacula using more invasive methods (Cowles, 1941; Sage, 2005; Burger et al., 2012) can directly change the behaviors of hibernating snakes or alter the physical conditions of the site.

During hibernation, snakes survive by freeze avoidance or supercooling and exhibit variable tolerances to water inundation and hypoxia. Snakes generally avoid freezing by retreating into a thermally buffered subterranean space in the autumn and remaining there until spring (Costanzo and Claussen, 1988; Storey and Storey, 1992; Gienger and Beck, 2011). Snakes are freeze-intolerant (Costanzo and Claussen, 1988; Storey and Storey, 1996; Costanzo and Lee, 2013), even species with very northerly distributions (Churchill and Storey, 1992; Lee and Constanzo, 1998; Andersson and Johansson, 2001; Costanzo and Lee, 2013). Supercooling (Storey and Storey, 1996) is only a short-term strategy for snakes at high latitude (Storey and Storey, 1992; Andersson and Johansson, 2001). There is some indirect evidence that hibernating snakes may drown if the hibernaculum floods (Viitanen, 1967; Prestt, 1971; Shine and Mason, 2004; Harvey et al., 2014), although access to water during overwintering can be beneficial to overall survival (Costanzo, 1989a,b; Prior, 1997; Costanzo et al., 2001). Snakes may survive flooding events in the field (May et al., 1996; Seigel et al., 1998) and in the laboratory (Constanzo, 1989a; Todd et al., 2009). However, the ability to remain completely submerged for the entire duration of hibernation seems unlikely and may vary among species and within species depending on hibernacula site quality (Gregory, 1982; Ultsch, 1989; Rollinson et al., 2008).

The survival of overwintering ectothermic vertebrates depends on the reduction in metabolic rates when exposed to cold environmental temperatures (reduced O₂ demand) and the ability to obtain oxygen via cutaneous respiration (Costanzo, 1989a; Tattersall and Boutilier, 1997; Schulte, 2015). Snakes in hibernacula may experience poor air quality when flooding, snow, or ice block atmospheric oxygen exchange (Boutilier, 1990), particularly when the surrounding organic soils produce toxic biogases (Moore and Knowles, 1989; Charman et al., 1994). Measurement of hibernaculum air quality is challenging, although assessment of groundwater dissolved oxygen (GWDO) could provide an indirect measure of aerobic quality changes within a hibernaculum.

We studied Massasauga (Sistrurus catenatus) hibernation habitat using a subterranean grid of wells and frost tubes at two sample areas (Mined and Not Mined), a system that allowed us to avoid disturbing the ecological functions of the hibernation site or the hibernating snakes themselves. A 4-yr radiotelemetry study (i.e., 2000-2004) averaging 5.2 Massasaugas per year with 5 individuals followed over multiple years (Yagi and Tervo, 2005) demonstrated that (1) individuals moved between Mined and Not-Mined habitats during the active season, (2) individuals exhibited site fidelity to either Mined or Not-Mined hibernation sites across years, and (3) not all hibernacula could be identified. To survive, we assumed that snakes would move within a vertical subterranean space that is aerobic and flood and frost free (Fig. 1), and we named this space the “life zone” (LZ). We measured groundwater level (GWL), frost depth (FD), groundwater temperature (T_gw), and dissolved oxygen (GWDO) at each site and hypothesized that (1) if an LZ requires a frost- and flood-free subterranean space, then we expect the LZ size to be affected by environmental variables, such as precipitation and temperature that cause flooding and freezing; (2) if the thermal and aerobic quality of the subterranean space is affected by the size of the LZ, then we predict an increasing quality with LZ size; and (3) if successful snake hibernation requires the presence of an LZ, then we expect to find increased snake overwintering survival with larger LZ size and higher quality. Additionally, we propose a standardized method for measuring, delineating, and communicating the characteristics of the LZ within terrestrial and semi terrestrial habitats using a long-term case study of an endangered Massasauga (Sistrurus catenatus) population.

Fig. 1.

Life-zone (LZ) model hypothesis for terrestrial and semiterrestrial hibernating Massasaugas. Life zone is the vertical subterranean space that remains aerobic, not frozen, and flood free throughout hibernation.

Materials and Methods

Study Site.—The study site is a partially drained and strip-mined peatland wetland ecosystem within the Great Lakes lowland region near Port Colborne Ontario, Canada (42°52′50′'N, 79°15′00′'W, datum WGS 84). The wetland is approximately 1,700 ha containing naturalized and remnant swamp, marsh, and bog vegetation communities within an organic soil basin underlain by thick clay soils (Browning, 2015). Constructed watercourses (i.e., drains) intersect both the organic and underlying clay layer; thus, shortening the hydroperiod (Browning, 2015; Fig. 2). When peat mining ceased in the 1990s, increases in internal drain water levels occurred because of intrinsic and extrinsic factors including vegetation growth, sediment accumulations, woody debris, constructed peat dams, Beaver (Castor canadensis) dams (Yagi and Litzgus, 2012; Browning, 2015), and stochastic wet weather events (Environment and Climate Change Canada, 2016). The wetland is isolated by surrounding agricultural, industrial, and rural land use and is more than 200 km from other Massasauga populations.

Fig. 2.

The Massasauga population study area is a partially strip-mined peatland located in southern Ontario, Canada. (A) shows the amount of surface water present during dry weather cycles (1998–2005; 2011–2014). (B) shows the amount of surface water present during a wet weather cycle (2006–2011). The first known flood event occurred from 10 October 2006 to 13 October 2006. Hibernation study areas are data sensitive.

Measurement of the LZ.—We initiated our LZ study near the end of the initial radiotelemetry study (Yagi and Tervo, 2005) sampling between the winter of 2003–2004 through the winter of 2013–2014. Mined and Not-Mined LZ study areas were chosen based upon hibernation site fidelity data. More than 1 km of forested peatland separates the LZ study areas. The study areas are also within different hydrologic subcatchments (Browning, 2015) but otherwise have a similar vegetation community (Betula pendula, Aronia melanocarpa, Vaccinium sp., Rubus sp.) with holes at the surface. The Mined site is low-lying with exposed organic soil surface flattened and compacted from past peat mining activities. The Not-Mined site is higher in elevation with more microtopography.

Groundwater level (GWL [cm]), FD (cm), groundwater temperature (T_gw [°C]), GWDO (mg/L), and snow height (cm) were measured manually in a well-field grid within the Mined and Not-Mined hibernation study areas for 11 winters. Winters were defined annually from 1 November to 30 April. Detailed methods for constructing and sampling the LZ well-field grid are in Appendix 1. The LZ (cm) calculation equals the difference between the top-of-well-pipe water-level (T) and the height of well above ground (H) plus FD (LZ = T – [H + FD]); or simply LZ = GWL – FD (Fig. 3). Observations of surface flooding or ice formation in hibernation areas constituted a zero LZ. We calculated the LZ groundwater thermal buffering function (T_buff °C) for each study area as the difference between T_gw and daily mean air temperature (T_AIR), or T_buff = T_gw – T_AIR (Environment and Climate Change Canada, 2016).

Fig. 3.

Hibernation habitat study methods developed for a partially mined peatland ecosystem inhabited by a Massasauga population located in southern Ontario, Canada. Life-zone methods include a grid of groundwater wells paired with a frost tube in each study area (Mined and Not Mined). Life zone (LZ) equals the groundwater level (GWL) minus frost depth (FD). All measurements are relative to the ground surface.

An important hydrological event occurred 3 yr into the LZ study when the Mined study area (0.8 ha) and approximately 317 ha within the central peat barrens flooded for the first time (Fig. 2; Yagi and Litzgus, 2012). The first flood event occurred during a stochastic storm from 10 October 2006 to 13 October 2006 that included mixtures of snow, rain, and freezing temperatures (Environment and Climate Change Canada, 2016). The flooded area accounted for over 30% of the central peat barrens, and the Mined study area continued to experience flood events over the next five winters (i.e., winter of 2006–2007 to winter of 2010–2011). The flood events were followed by a dry weather cycle with a severe drought and central peatland wildfires in 2012 (Fig. 2).

Because we were unable to maintain a strict weekly sampling schedule over the 11-yr study period, and because field equipment sometimes failed, temporal gaps in our data set exist. For missing GWL, we used the average between subsequent measurements and assumed the levels did not change drastically between visits. We used the average daily FD value for the study area to fill weekly gaps to complete the LZ calculation. Although our approach had a smoothing effect on the data, we were able to retain the spatial context. The data correction occurred 146 times out of 3170 (4.6%). We omitted blank GWDO and T_gw data from further analysis.

Massasauga Population.—We used annual (April–November) mark–recapture (M-R) data from 1998 to 2016 to monitor population trends (Ontario Ministry of Natural Resources and Forestry, unpubl. data). All captured adults were tagged subcutaneously with a passive integrated transponder (PIT; Parent, 1997). We identified neonatal snakes by their head and mottle pattern. We estimated age from several biological parameters such as snout–vent length (SVL), weight, reproductive maturity, and the number of rattles within a complete rattle sequence. We assumed that two rattles form per year for this population, which is consistent with local data and other rattlesnake populations at similar latitudes (Fitch, 1985; Macartney et al., 1990; Aldridge and Brown, 1995). All encountered snakes were aged in this manner to confirm the presence in the population and allow us to discern a preflood group from postflood recruits. The number of snakes considered to be adults (N_adults) was aged 3 yr and older.

We did not use traditional calculations for population estimates because of the cryptic nature of Massasaugas, our low annual recapture rates, and uncertainty in annual survival because of the presence of impaired habitat and environmental stochasticity. We instead built a time series of known individuals (N_indiv) for the overall sampling area and LZ study period. All encountered snakes per year plus those presumed present because of their age estimates from future encounters (N_indiv = new captures [C] + recaptures [R] + known undetected [U]). We calculated the yearly rate of detection (D) as (C + R)/N_indiv. Snakes not recaptured up to 10 yr after the last encounter were hence removed from the sample population. We calculated the proportion of adults in the population (P_adults) as N_adults/N_indiv.

Statistical Analyses.—We used the program R to complete all statistical analyses and graphing (version 3.5.0; R Core Team, 2018). We tested parameters and model residuals for normality using a Shapiro-Wilk test and visually using qq-plots. We used nonparametric tests (i.e., Kruskal-Wallis and Dunn's test) when transformations failed to correct for violations of parametric statistical assumptions.

Weather Data Analysis.—We downloaded 16 winters (1998–2014) of historical weather data from the nearest weather station with vetted data (Port Colborne or Hamilton, Ontario). This period encompassed the broader 16-yr Massasauga population study, including the 11-yr LZ study. Thus, it allowed us to contextualize the LZ study within recent winter climate trends concerning precipitation and temperature. We partitioned the monthly weather data to match the LZ study period and extracted four weather variables for analysis: precipitation (snow, rain, total) and the number of freezing days (i.e., mean daily temperature < 0°C; Environment and Climate Change Canada, 2016). We conducted a principal component analysis (PCA) to reduce the winter weather records into useful categories for analysis. We employed a K-means cluster analysis with the K-means function to assess the sums of squares based on the different number of clusters. By plotting the total within-groups sums of squares against the number of clusters, the sudden change in slope depicted an optimal value of three winter types (Everitt et al., 2011). We partitioned the three winter types into groups defined by temperature and precipitation trends: cold–dry (winter type 1), mild–wet (winter type 2), and cold–wet (winter type 3). We achieved a data reduction of the same weather data using PCA with the FactoMineR package (Lê et al., 2008). All variables were centered and scaled to unit variance before the PCA. The predominant axis (PC1) described a measure of rainfall (from dry to wet), and the second axis (PC2) combined the number of days below 0°C (from mild to cold) and total precipitation as snow.

Life Zone Analysis.—The LZ is the space between the FD and GWL. Therefore, according to our hypothesis, when GWL exceeds the ground surface or frost reaches the GWL, LZ size equals zero. However, we used negative LZ values to avoid zero truncation issues during the model analysis (Zuur et al., 2009). Fixed and random effects (Table 1) were analyzed with respect to the response variables LZ size and natural log-transformed GWDO and T_buff using linear mixed effect models on the full data set with the lmer function from the lmerTest package. Adjusted degrees of freedom were estimated using the Satterthwaite method. A model selection analysis using the AIC_c function from the AICcmodavg package was used to assess the fixed effects that best predict LZ size, ln GWDO, and T_buff (Akaike, 1987). A Type III ANOVA was used to analyze the significance between fixed effects within the Akaike top model.

Table 1.

Summary of parameters used in the model selection analysis for each habitat and encounter response variable from an 11-yr hibernation habitat study on a Massasauga population located in southern Ontario, Canada. Linear mixed-effects models are (1) Life-zone size (LZ [cm]), (2) LZ thermal buffering function (T_buff), (3) groundwater dissolved oxygen (GWDO [mg/L]), (4) a generalized linear model compares the previous winter's fixed effects against the following active-season Massasauga encounters (N_indiv).

Massasauga Survival Analysis.—We used a linear model to assess differences in detection by flood groupings. An annual flood survival estimate was calculated for the preflood group accounting for known deaths. We used a generalized linear model with a Poisson distribution to compare the previous winter's fixed effects (LZ size, T_buff, GWDO, T_gw, snow height) against the following active season N_indiv (Table 1). Akaike model selection was used to determine the top model. A Type II analysis of deviance was used to assess the significance of the fixed-effects within the Akaike top model, using the ANOVA function from the R car package.

Results

Winter Weather Trends.—The PCA combined local winter environmental conditions into PC1 (rainfall) and PC2 (snow and number of freezing days), which clustered into three winter types (Fig. 4; Table 2). Winter type 1 (cold–dry) occurred four times but only once during the LZ study. Winter type 2 (mild–wet) occurred four times during the LZ study. Winter type 3 (cold-wet) occurred six times (Fig. 4).

Fig. 4.

PCA plot depicting the winter-type categories developed for winter environmental data for a Massasauga population located in southern Ontario, Canada. We analyzed local annual winter (November–April) environmental data: total precipitation, rainfall, snow, and the number of days mean temperature is below zero (Environment and Climate Change Canada, 2016).

Table 2.

Summary of PCA eigenvalues for winter environmental data from a Massasauga population in southern Ontario, Canada. We downloaded historical weather data for the study area from the Environment and Climate Change Canada weather data website (1998–2014). Variables of interest include total rain, total snow, total precipitation, and the number of days when the mean temperature was below zero.

LZ Size Dynamics.—The Not-Mined area has significantly deeper organic soils (mean = 167.9 cm, SD = 53.3) vs. Mined (mean = 71.3 cm, SD =20.7; t₍₃₀₎ = 8.49, P < 0.001). Maximum recorded frost depth (29.8 cm) was in Not-Mined habitat during winter 2007 after the first flood event. The Mined LZ size was significantly larger before than after flooding (P < 0.001; Table 3). Mined LZ size was reduced to zero periodically (312/658 measurements = 47% zeros) during flood events (Fig. 5). Not-Mined LZ size was not significantly different throughout the 11-yr study (P = 0.49; Table 3; Fig. 5). The model that best explained the variance in LZ size included the fixed effects of site soil depth (cm) and flood condition (Table 4). Type III ANOVA indicated that site soil depth and flood condition were significant effects on LZ size (Table 5).

Table 3.

Life-zone summary measurements mean ± SD from Mined and Not-Mined study areas before and after first flood event which initiated a 5-yr stochastic flooding period followed by receding floodwaters (asterisk denotes significant P value).

Fig. 5.

Massasauga hibernation habitat weekly life-zone size (LZ cm) by study site (Mined and Not Mined), over time (November 2003–April 2014). The study site is a partially strip-mined peatland located in southern Ontario, Canada. Life-zone size was zero for at least 1 wk during winters (2006–2011) because of surface flooding and freezing during a wet weather cycle. Error shown shaded in gray is ±95% confidence interval (CI).

Table 4.

Akaike top five linear mixed effect models and the null model results for an 11-yr Massasauga hibernation habitat study located in southern Ontario, Canada. The models are (A) Life-zone size (LZ [cm]), (B): thermal buffering function (T_buff), and (C) groundwater dissolved oxygen (ln GWDO). Fixed effects include winter type (1–3) or PC1, or PC2, site (Mined or Not Mined) or site soil depth (cm), flood (before and after) or weather cycles (dry and wet), and LZ size (cm). For all analyses, random effects were Well ID, date, and year (null model). The top models are marked by an asterisk.

Table 5.

The top-ranked linear mixed-effects model ANOVA results for the 11-yr Massasauga hibernation habitat study located in southern Ontario, Canada. The hibernation habitat study examined the following environmental effects: (A) life-zone size (LZ cm), (B) LZ thermal buffering function (T_buff), and (C) LZ aerobic quality function (GWDO). Random effects were well ID, date, and year (* denotes significant P value).

Groundwater Temperature.—Mined T_gw declined significantly following flood events (P < 0.001), whereas Not-Mined T_gw did not (P = 0.08; Table 3). T_gw was not measurable when a well was flooded or frozen (361 out of 3,170 sampling events, 11.4%), causing data gaps and imbalances in the Mined area data set during the flood maxima, winter 2008. There were also 75 out of 3,170 (2.3%) sampling events that were too dry to measure groundwater attributes. Mined T_buff was not significantly different due to flooding (P = 0.11; Table 3); and Not-Mined T_buff was significantly greater after flooding (P < 0.001; Table 3). The top model for T_buff included the significant fixed effects of LZ size, snow height, site soil-depth, flood-condition, and PC1 (Tables 4 and 5).

Groundwater Dissolved Oxygen.—Mined GWDO (mg/L) was significantly reduced when flooded (P < 0.001) and Not-Mined GWDO significantly increased (P < 0.001; Table 3). Mined ln GWDO was significantly and positively correlated with mean LZ size before flooding (F_1,34 = 41.25, R² = 0.55, P < 0.001) but not after (F_1,41 = 0.48, R² = 0.01, P = 0.49; Fig. 6). Although the Not-Mined area did not flood, it did undergo the same wet and dry weather cycles. The relationship between Not-Mined mean GWDO and mean LZ size was significant before (F_1,32 = 2.6, R² = 0.14, P = 0.03) but not after (F_1,51 < 0.001, R² < 0.001, P = 0.99) the first flood event (Fig. 6). The top model that best explained the variance in ln GWDO included the fixed effects LZ size, PC1, flood condition, and snow height (Table 4). Type III ANOVA for the top model showed that LZ and PC1 were significant factors, whereas snow height and flood condition were not (Table 5).

Fig. 6.

Massasauga hibernation habitat weekly groundwater dissolved oxygen (mean GWDO mg/L) by weekly life-zone size (mean LZ cm). The study area is a partially mined peatland located in southern Ontario, Canada. Mined and Not-Mined areas are grouped by before (2003–2005) and after (2006–2014) first flood event. Error shown shaded in gray is ±95% CI.

Massasauga Population Dynamics.—Mean snake detection (D) increased in the study area from before (mean = 0.30, SD = 0.11) to during (mean = 0.52, SD = 0.37) and decreased after (mean = 0.38, SD = 0.19) the flood events. There was no significant difference in D among these flood groupings (F_2,8 = 0.68, P = 0.53). Therefore, we used before and after the first flood event groups for subsequent analyses. Total search effort during the 11-yr study averaged 735 person hours per year, SD = 354. The lowest N_indiv occurred in 2008, and the preflood group low occurred in 2010 (Fig. 7). The preflood group encounters declined by 33% after the first flooded winter, followed by 76% after the second, and 90% by the third winter (Fig. 7). The mean age of N_indiv in 2004 was 2.5 yr ± 0.50 SE (P_adults = 0.45), increased to 4.1 yr ± 0.46 S.E. in 2007 (P_adults = 0.79), and declined to 1 yr ± 0.6 S.E. by 2010 (P_adults = 0.11; Appendix 2). The Akaike top model that predicted N_indiv included LZ size over the fixed effects T_buff, snow height, T_gw, and GWDO (Table 6). The N_indiv and previous winter's mean LZ was a significant positive relationship with an incidence response ratio of 0.06 (χ² = 28.8, df =1, r² = 0.34, P < 0.001; Fig. 8).

Fig. 7.

Massasauga encounters (N_indiv) from 2004–2014 are calculated from the long-term (1998–2016) encounter data set for an isolated population located in southern Ontario, Canada. N_indiv includes all encountered adult and juvenile Massasaugas plus those presumed present due to their age estimates from future encounters. The preflood group encounter trend is indicated to discern postflood recruitment.

Table 6.

Akaike top five generalized linear mixed effect models and the null model results for an 11-yr Massasauga hibernation habitat and population study located in southern Ontario, Canada. Massasauga encounters (N_indiv) were set as a function of the fixed-effects mean; Life-zone size (LZ cm), snow height (cm), LZ temperature buffering function (T_buff), and LZ aerobic quality (GWDO). N_indiv includes all encountered snakes per year plus those presumed present because of their age estimates from future encounters.

Fig. 8.

Generalized linear regression models are provided for Massasauga encounters (N_indiv) with respect to the independent measures of the previous winter's mean annual, (A) groundwater thermal buffer index (T_buff = (T_gw – Mean daily T_Air ), (B) groundwater dissolved oxygen (GWDO mg/L), (C) snow height (cm), and (D) the Akaike top model is Life-zone size (LZ cm). Error shown shaded gray is ± SE. Encounter data are from (1998–2016), for an isolated Massasauga population located in southern Ontario, Canada.

Discussion

In our study, anthropogenic habitats and natural habitats show important differences in life-zone functions during times of environmental stochasticity. Important factors affecting LZ size in our model analysis were site, soil depth, snow height, flood condition, and winter type. The significant relationship between LZ size and flood condition supported our first hypothesis that environmental factors, such as precipitation and temperature, affect the size of the LZ. Thermal buffering and GWDO quality attributes significantly improved in natural areas that maintained an LZ compared to peat-mined habitats that did not. The significant relationship between LZ size, T_buff, and GWDO supports our second hypothesis that LZ size affects thermal and aerobic quality. Finally, the significant relationship between LZ size and Massasauga encounters indirectly supports our third hypothesis that Massasauga survival increases with LZ size. Because natural habitats only maintained an LZ, they likely provided areas of refugia for the Massasauga population.

In our LZ model, we considered that snakes might die from either drowning, freezing, or asphyxiation. However, only the top model showed a positive correlation between N_indiv and LZ size. In altered ecosystems, a large LZ may provide the capacity needed to maintain survival during times of environmental stochasticity. In natural bog ecosystems, a shallow LZ may also provide thermal stability because there is an elevated and stable groundwater table (Smolarz et al., 2018). At our site, we have a degraded ombrotrophic bog ecosystem without stable hydrology or stable Massasauga population. Therefore, relationships between population trends and LZ functions may be confounded by these complexities. In the case of GWDO, we did not directly measure the air oxygen content in occupied hibernacula. Yet, over half the measurements of GWDO were hypoxic, which indicates that cutaneous respiration, as a flood survival strategy, may not be possible here (Costanzo and Lee, 1995; Tattersall and Boutilier, 1997; Jackson, 2007). Therefore, maintenance of an aerobic space within the hibernation site for oxidative metabolism is an important consideration in overwintering survival. Because subterranean air spaces naturally have low atmospheric oxygen exchange, winter types vary, and flooding reduces soil oxygen content, additional monitoring is needed to establish the relationship between LZ size, survival, GWDO, snow height, and air space O₂ (Boutilier, 1990; Cavallaro and Hoback, 2014). Survival experiments that control snake overwintering locations, coupled with LZ measures, would provide direct support of our hypothesis that the presence of an LZ supports overwinter survival.

Our population trends are based upon Massasauga encounters that are backcast using age and refined annually with new M-R and age estimate data. Therefore, N_indiv is more reliable in the past and underestimated in the most recent years. Our M-R data are affected by the challenges of finding a cryptic species. Yet our snake detection rate doubled during the first three flood years, then dropped to zero in the fourth flood year (i.e., 2010). Zero captures in 2010 are an important finding, considering there was a spatial search effort bias toward the central peat barrens as we were trying to recapture flood survivors. After 2010 we evened out the search effort and captured Massasaugas in the Not-Mined habitats that were previously undetected. Declining encounters and the significant correlation with LZ size provide indirect evidence of reduced survival in the wetland during this period of environmental stochasticity. We also presented evidence that the mean age of N_indiv increased by a factor of 1.6 from 2004 to 2007, then declined by a factor of 4.1 by 2010. The increase in the mean age estimate suggests (1) low recruitment or reproductive output, (2) that adults had a higher initial survival rate than their offspring, or (3) a combination of factors. The subsequent decline in mean age by 2010 is likely related to successful overwintering in the younger age classes. Because postflood encounters were almost exclusively in Not-Mined habitats that maintained an LZ, females giving birth in higher elevations may have aided neonatal dispersal and survival. Adults likely moved to higher-elevation areas and overwintered successfully, emigrated beyond the survey areas, or did not survive hibernation in low-lying areas because of site fidelity. This study presents indirect evidence of the latter. With recent flooding events and the challenges in monitoring a recovering cryptic species, the relationship between LZ functions and survival requires further investigation. We should avoid overinterpreting these results because our summarized data may mask the true relationship between encounters and hibernation habitat quality.

Across the Massasauga range, hibernation sites are often, but not always, found in habitats with an elevated groundwater table (Johnson, 1995; Parent, 1997; Harvey and Weatherhead, 2006b). Habitats that maintain an LZ may offer increased thermal stability and population viability (Pomara et al., 2014; Smolarz et al., 2018). The presence of a consistent snow layer may also provide important thermal buffering, especially when considering the presence of a thermally stable groundwater table and the thickness of the LZ. In addition to the increased frequency of stochastic events that reduce LZ size, the presence of a snow layer at our study site in southern Ontario is inconsistent between years. Increasingly episodic polar vortices (Charlton and Polvani, 2007) may also increase freeze risk in areas, especially when there is a lack of snow cover and a small LZ. Modeling LZ size and quality functions within various climate-change scenarios under ecologically relevant conditions are important next steps in this research (Pomara et al., 2014). Expanding LZ monitoring into other habitats and soil types are warranted where there is a flood or freeze concern.

Furthermore, demonstrating causality by multiple potential factors is challenging, especially when subsurface mortality and unknown hibernacula locations prevent the confirmation of survival. The continual presence of potentially poor habitat that is periodically used by Massasaugas during dry weather cycles may indicate the presence of an ecological trap (Battin, 2004). From a conservation perspective, mortality events caused by an ecological trap may increase the extinction risk of a population that may already be compromised by small size, lack of rescue by isolation, genetic drift, and uncertainty surrounding adaptation to environmental stochasticity (Battin, 2004; Bradke et al., 2018). Future studies on this topic should provide opportunities to study the effects of climate change, habitat quality, winter severity, and the relevance to reptile hibernation success in northern climates.