Abstract

Temperate snakes occupy overwintering sites for most of their annual life cycle. Microhabitat characteristics of the hibernaculum are largely undescribed, yet are paramount in ensuring snake overwintering survival. We hypothesized that snakes survive hibernation within a vertical subterranean space that we termed a “life zone” (LZ), that is aerobic and flood and frost free throughout winter. We studied an isolated, endangered population of Massasaugas (Sistrurus catenatus) inhabiting an anthropogenically altered peatland and monitored the subterranean habitat during a period of environmental stochasticity. Initial radio telemetry confirmed that snakes moved between altered and natural habitats during the active season and showed hibernation-site fidelity to either habitat. We used a grid of groundwater wells and frost tubes installed in each hibernation area to measure LZ characteristics over 11 consecutive winters. The LZ within the impacted area was periodically reduced to zero during a flood–freeze cycle, but the LZ in the natural area was maintained. Model selection analysis revealed that soil depth and flood status best predicted LZ size. Thermal buffering and groundwater dissolved oxygen increased with LZ size, and annual Massasauga encounters were significantly correlated with LZ size. This analysis suggests a population decline occurred when LZ size was reduced by flooding. Our data give support to the importance and maintenance of an LZ for successful snake hibernation. Our methods apply to subterranean hibernation habitats that are at risk of environmental stochasticity, causing flooding, freezing, or hypoxia.

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 O2 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 (TGW), 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.

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.

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 (TGW [°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 (TBUFF °C) for each study area as the difference between TGW and daily mean air temperature (TAIR), or TBUFF = TGW − TAIR (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.

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 TGW 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 (Nadults) 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 (Nindiv) 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 (Nindiv = new captures [C] + recaptures [R] + known undetected [U]). We calculated the yearly rate of detection (D) as (C + R)/Nindiv. 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 (Padults) as Nadults/Nindiv.

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 TBUFF 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 AICc function from the AICcmodavg package was used to assess the fixed effects that best predict LZ size, ln GWDO, and TBUFF (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 (TBUFF), (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 (Nindiv).

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 (TBUFF), (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 (Nindiv).
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 (TBUFF), (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 (Nindiv).

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, TBUFF, GWDO, TGW, snow height) against the following active season Nindiv (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).

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.

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.
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(30) = 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).

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

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 (TBUFF), 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.

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 (TBUFF), 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.
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 (TBUFF), 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 (TBUFF), and (C) LZ aerobic quality function (GWDO). Random effects were well ID, date, and year (* denotes significant P value).

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 (TBUFF), and (C) LZ aerobic quality function (GWDO). Random effects were well ID, date, and year (* denotes significant P value).
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 (TBUFF), and (C) LZ aerobic quality function (GWDO). Random effects were well ID, date, and year (* denotes significant P value).

Groundwater Temperature.—

Mined TGW declined significantly following flood events (P < 0.001), whereas Not-Mined TGW did not (P = 0.08; Table 3). TGW 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 TBUFF was not significantly different due to flooding (P = 0.11; Table 3); and Not-Mined TBUFF was significantly greater after flooding (P < 0.001; Table 3). The top model for TBUFF 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 (F1,34 = 41.25, R2 = 0.55, P < 0.001) but not after (F1,41 = 0.48, R2 = 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 (F1,32 = 2.6, R2 = 0.14, P = 0.03) but not after (F1,51 < 0.001, R2 < 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.

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 (F2,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 Nindiv 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 Nindiv in 2004 was 2.5 yr ± 0.50 SE (Padults = 0.45), increased to 4.1 yr ± 0.46 S.E. in 2007 (Padults = 0.79), and declined to 1 yr ± 0.6 S.E. by 2010 (Padults = 0.11; Appendix 2). The Akaike top model that predicted Nindiv included LZ size over the fixed effects TBUFF, snow height, TGW, and GWDO (Table 6). The Nindiv and previous winter's mean LZ was a significant positive relationship with an incidence response ratio of 0.06 (χ2 = 28.8, df =1, r2 = 0.34, P < 0.001; Fig. 8).

Fig. 7

Massasauga encounters (Nindiv) from 2004–2014 are calculated from the long-term (1998–2016) encounter data set for an isolated population located in southern Ontario, Canada. Nindiv 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.

Fig. 7

Massasauga encounters (Nindiv) from 2004–2014 are calculated from the long-term (1998–2016) encounter data set for an isolated population located in southern Ontario, Canada. Nindiv 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 (Nindiv) were set as a function of the fixed-effects mean; Life-zone size (LZ cm), snow height (cm), LZ temperature buffering function (TBUFF), and LZ aerobic quality (GWDO). Nindiv includes all encountered snakes per year plus those presumed present because of their age estimates from future encounters.

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 (Nindiv) were set as a function of the fixed-effects mean; Life-zone size (LZ cm), snow height (cm), LZ temperature buffering function (TBUFF), and LZ aerobic quality (GWDO). Nindiv includes all encountered snakes per year plus those presumed present because of their age estimates from future encounters.
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 (Nindiv) were set as a function of the fixed-effects mean; Life-zone size (LZ cm), snow height (cm), LZ temperature buffering function (TBUFF), and LZ aerobic quality (GWDO). Nindiv 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 (Nindiv) with respect to the independent measures of the previous winter's mean annual, (A) groundwater thermal buffer index (TBUFF = (TGW – Mean daily TAir ), (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.

Fig. 8

Generalized linear regression models are provided for Massasauga encounters (Nindiv) with respect to the independent measures of the previous winter's mean annual, (A) groundwater thermal buffer index (TBUFF = (TGW – Mean daily TAir ), (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, TBUFF, 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 Nindiv 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 O2 (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, Nindiv 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 Nindiv 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.

Acknowledgments

All procedures involving animals followed the Ontario Ministry of Natural Resources and Forestry (OMNRF) regulations, and Animal Care permit 98–55 to 16–55. We have withheld locations of hibernation study areas because of data sensitivity. Environment Canada Habitat Stewardship Fund, OMNRF, and the World Wildlife Fund provided funding support. We thank our partners, Land Care Niagara, Niagara Peninsula Conservation Authority, National Massasauga Recovery Team, Toronto Zoo, B. Johnson, J. Litzgus, Vittie family, Yagi family, J. Durst, M. Browning, K. Frohlich, D. Mills, A. Parks, D. Denyes, M. Karam, and C. Blott. This article is dedicated to the memory of R. Tervo (d. 8 November 2008) and R. J. Planck (d. 29 October 2015).

Literature Cited

Literature Cited
Agugliaro,
J.
2011
.
Ecophysiology of hibernation in timber rattlesnakes (Crotalus horridus)
.
Ph.D. diss.
,
University of Arkansas
,
Arkansas, USA
.
Akaike,
H.
1987
.
Factor analysis and AIC
.
Psychometrika
52
:
317
332
.
Aldridge,
R. D.,
and
W. S.
Brown.
1995
.
Male reproductive cycle, age at maturity, and cost of reproduction in the timber rattlesnake (Crotalus horridus)
.
Journal of Herpetology
29
:
399
407
.
Andersson,
S.,
and
L.
Johansson.
2001
.
Cold hardiness in the boreal adder, Vipera berus
.
Cryo-Letters
3
:
151
156
.
Battin,
J.
2004
.
When good animals love bad habitats: ecological traps and the conservation of animal populations
.
Conservation Biology
18
:
1482
1491
.
Boutilier,
R. G.
1990
.
Respiratory gas tensions in the environment
.
Advances in Comparative and Environmental Physiology
6
:
1
13
.
Bradke,
D. R.,
E. T.
Hileman,
J. F.
Bartman,
L. J.
Faust,
R. B.
King,
N.
Kudla,
and
J. A.
Moore.
2018
.
Implications of small population size in a threatened pitviper species
.
Journal of Herpetology
52
:
387
397
.
Browning,
M. H. R.
2015
.
The dynamics and mechanisms of community assembly in a mined Carolinian peatland
.
Ph.D. diss.
,
Trent University
,
Peterborough, Canada
.
Burger,
J.,
R. T.
Zappalorti,
M.
Gochfield,
E.
DeVito,
D.
Schnieder,
and
C.
Jeitner.
2012
.
Long-term use of hibernacula by northern pinesnakes (Pituophis melanoleucus)
.
Journal of Herpetology
60
:
596
601
.
Cavallaro,
M. C.,
and
W. W.
Hoback.
2014
.
Hypoxia tolerance of larvae and pupae of the semi-terrestrial caddisfly (Trichoptera: Limnephilidae)
.
Ecology and Population Biology
107
:
1081
1085
.
Charlton,
A. J.,
and
L. M.
Polvani.
2007
.
A new look at stratospheric sudden warmings. Part I: climatology and modeling benchmarks
.
Journal of Climate
20
:
449
469
.
Charman,
D. J.,
R.
Aravena,
and
B. G.
Warner.
1994
.
Carbon dynamics in a forested peatland in north-eastern Ontario, Canada
.
Journal of Ecology
82
:
55
62
.
Churchill,
T. A.,
and
K. B.
Storey.
1992
.
Freezing survival of the gartersnake Thamnophis sirtalis parietalis
.
Canadian Journal of Zoology
70
:
99
105
.
Costanzo,
J. P.
1989
a
.
A physiological basis for prolonged submergence in hibernating gartersnakes Thamnophis sirtalis: evidence for an energy-sparing adaptation
.
Physiology of Zoology
62
:
580
592
.
Costanzo,
J. P.
1989
b
.
Effects of humidity, temperature, and submergence behaviour on survivorship and energy use in hibernating gartersnakes, Thamnophis sirtalis
.
Canadian Journal of Zoology
67
:
2486
2492
.
Costanzo,
J. P.,
and
D. L.
Claussen.
1988
.
Natural freeze tolerance in a reptile
.
Cryo-Letters
9
:
380
385
.
Costanzo,
J. P.,
and
R. E.
Lee
Jr.
1995
.
Supercooling and ice nucleation in vertebrate ectotherms
.
Pp
.
221
237
in
R. E.
Lee,
G. J.
Warren,
and
L.V.
Gusta
(
eds.
),
Biological Ice Nucleation and Its Applications
.
APS Press
,
USA
.
Costanzo,
J. P.,
J. D.
Litzgus,
J. B.
Iverson,
and
R. E.
Lee
Jr.
2001
.
Cold-hardiness and evaporative water loss in hatchling turtles
.
Physiological and Biochemical Zoology
74
:
510
519
.
Costanzo,
J. P.,
and
R. E.
Lee
Jr.
2013
.
Commentary: avoidance and tolerance of freezing in ectothermic vertebrates
.
The Journal of Experimental Biology
216
:
1961
1967
.
Cowles,
R. B.
1941
.
Observations on the winter activities of desert reptiles
.
Ecology
22
:
125
140
.
Environment and Climate Change Canada
.
2016
.
Historical climate data; 1998 to 2014, Port Colborne and Hamilton, Canada
.
Everitt,
B. S.,
S.
Landau,
M.
Leese,
and
D.
Stahl.
2011
.
Hierarchical clustering
.
Pp
.
71
110
in
B. S.
Everitt,
S.
Landau,
M.
Leese,
D.
Stahl
(
eds.
),
Cluster Analysis. 5th ed
.
John Wiley & Sons
,
Inc., USA
.
Fitch,
H. S.
1985
.
Observations on rattle size and demography of prairie rattlesnakes (Crotalus viridis) and timber rattlesnakes (Crotalus horridus) in Kansas
.
Museum of Natural History, University of Kansas
118
:
1
11
.
Gienger,
C. M.,
and
D. D.
Beck.
2011
.
Northern Pacific rattlesnakes (Crotalus oreganus) use thermal and structural cues to choose overwintering hibernacula
.
Canadian Journal of Zoology
89
:
1084
1090
.
Gregory,
P. T.
1982
.
Reptilian hibernation
.
Pp
.
53
154
in
C.
Gans
and
F.H.
Pough
(
eds.
),
Biology of the Reptilia
.
Academic Press
,
USA
.
Harvey,
D. S.,
and
P. J.
Weatherhead.
2006
a
.
A test of the hierarchical model of habitat selection using eastern massasauga rattlesnakes (Sistrurus c. catenatus)
.
Biological Conservation
130
:
206
216
.
Harvey,
D. S.,
and
P. J.
Weatherhead.
2006
b
.
Hibernation site selection by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) near their northern range limit
.
Journal of Herpetology
40
:
66
73
.
Harvey,
D. S.,
A. M.
Lentini,
K. M.
Cedar,
and
P. J.
Weatherhead.
2014
.
Moving massasaugas: insight into rattlesnake relocation using Sistrurus c. catenatus
.
Herpetological Conservation and Biology
9
:
67
75
.
Huey,
R. B.
1982
.
Temperature, physiology, and the ecology of reptiles
.
Pp
.
25
91
in
C.
Gans,
and
F. H.
Pough
(
eds.
),
Biology of the Reptilia
.
Academic Press
,
USA
.
Jackson,
D. C.
2007
.
Temperature and hypoxia in ectothermic tetrapods
.
Journal of Thermal Biology
32
:
125
133
.
Johnson,
G.
1995
.
Spatial Ecology, Habitat Preference, and Habitat Management of the Eastern Massasauga, Sistrurus c. catenatus in a New York Weakly-Minerotrophic Peatland
.
Ph.D. diss.
,
State University of New York College of Environmental Science and Forestry
,
New York, USA
.
Jones,
P. C.,
R. B.
King,
R. L.
Bailey,
N. D.
Bieser,
K.
Bissell,
H.
Campa,
III,
T.
Crabill,
M. D.
Cross,
B. A.
Degregorio,
M. J.
Dreslik,
et al.
2012
.
Population ecology range-wide analysis of eastern massasauga survivorship
.
The Journal of Wildlife Management
76
:
1576
1586
.
Kingsbury,
B. A.,
and
C. J.
Coppola.
2000
.
Hibernacula of the copperbelly water snake (Nerodia erythrogaster neglecta) in southern Indiana and Kentucky
.
Journal of Herpetology
34
:
294
298
.
Lê,
S.,
J.
Josse,
and
F.
Husson.
2008
.
FactoMineR: an R package for multivariate analysis
.
Journal of Statistical Software
25
:
1
18
.
Lee,
R. E.,
Jr.,
and
J. P.
Costanzo.
1998
.
Biological ice nucleation and ice distribution in cold-hardy ectothermic animals
.
Annual Review of Physiology
60
:
55
72
.
Lentini,
A. M.,
G. J.
Crawshaw,
L. E.
Licht,
and
D. J.
McLelland.
2011
.
Pathologic and hematologic responses to surgically implanted transmitters in eastern massasauga rattlesnakes (Sistrurus catenatus catenatus)
.
Journal of Wildlife Diseases
47
:
107
125
.
Macartney,
J. M.,
K. W.
Larsen,
and
P. T.
Gregory.
1989
.
Body temperatures and movements of hibernating snakes (Crotalus and Thamnophis), and thermal gradients of natural hibernacula
.
Canadian Journal of Zoology
67
:
108
114
.
Macartney,
J. M.,
P. T.
Gregory,
and
M. B.
Charland.
1990
.
Growth and sexual maturity of the western rattlesnake, Crotalus viridis, in British Columbia
.
Copeia
1990
:
528
542
.
May,
P. G.,
T. M.
Farrell,
S. T.
Heulett,
M. A.
Pilgrim,
L. A.
Bishop,
D. J.
Spence,
A. M.
Rabatsky,
M. G.
Campbell,
A. D
Aycrigg,
and
W. E.
Richardson,
II.
1996
.
Seasonal abundance and activity of a rattlesnake (Sistrurus miliarius barbouri) in central Florida
.
Copeia
1996
:
389
401
.
Moore,
T. R.,
and
R.
Knowles.
1989
.
The influence of water table levels on methane and carbon dioxide emissions from peatland soils
.
Canadian Journal of Soil Science
69
:
33
38
.
Nordberg,
E. J.,
and
V. A.
Cobb.
2016
.
Midwinter emergence in hibernating timber rattlesnakes (Crotalus horridus)
.
Journal of Herpetology
50
:
203
208
.
Parent,
C. E.
1997
.
The Effects of Human Disturbance on Eastern Massasauga Rattlesnakes (Sistrurus catenatus catenatus) in Killbear Provincial Park
.
M.S. thesis
,
Carleton University
,
Canada
.
Patric,
J. H.,
and
B. D.
Fridley.
1969
.
A device for measuring soil frost
.
USDA Forest Service Research Note, U.S. Department of Agriculture
94
:
1
7
.
Pomara,
L. Y.,
O. E.
LeDee,
K. J.
Martin,
and
B.
Zuckerberg.
2014
.
Demographic consequences of climate change and land cover help explain a history of extirpations and range contraction in a declining snake species
.
Global Change Biology
20
:
2087
2099
.
Prestt,
I.
1971
.
An ecological study of the viper Vipera berus in southern Britain
.
Journal of Zoology
164
:
373
418
.
Prior,
K. A.
1997
.
Conservation Biology of Black Rat Snakes: Ecological, Demographic, and Genetic Approaches
.
Ph.D. diss.
,
Carleton University
,
Canada
.
R Core Team
.
2018
.
R: A language and environment for statistical computing
.
R Foundation for Statistical Computing
,
Vienna, Austria
.
Refsnider,
J. M.,
J.
Strickland,
and
F. J.
Janzen.
2012
.
Home range and site fidelity of imperiled ornate box turtles (Terrapene ornata) in northwestern Illinois
.
Chelonian Conservation and Biology
11
:
78
83
.
Rollinson,
N.,
G. J.
Tattersall,
and
R. J.
Brooks.
2008
.
Overwintering habitats of a northern population of painted turtles (Chrysemys picta): winter temperature selection and dissolved oxygen concentrations
.
Journal of Herpetology
42
:
312
321
.
Sage,
J. R.
2005
. Spatial ecology, habitat utilization, and hibernation ecology of the eastern massasauga (Sistrurus catenatus catenatus) in a disturbed landscape
.
M.S. thesis
,
Purdue University
,
Indiana, USA
.
Schulte,
P. M.
2015
.
The effects of temperature on aerobic metabolism: towards a mechanistic understanding of the responses of ectotherms to a changing environment
.
Journal of Experimental Biology
218
:
1856
1866
.
Seigel,
R. A.,
C. A.
Sheil,
and
J. S.
Doody.
1998
.
Changes in a population of an endangered rattlesnake Sistrurus catenatus following a severe flood
.
Biological Conservation
83
:
127
131
.
Shine,
R.,
and
R. T.
Mason.
2004
.
Patterns of mortality in a cold-climate population of gartersnakes (Thamnophis sirtalis parietalis)
.
Biological Conservation
120
:
201
210
.
Smolarz,
A. G.,
P. A.
Moore,
C. E.
Markle,
and
J. M.
Waddington.
2018
.
Identifying resilient eastern massasauga rattlesnake (Sistrurus catenatus) peatland hummock hibernacula
.
Canadian Journal of Zoology
,
96
:
1024
1031
.
Storey,
K. B.,
and
J. M.
Storey.
1992
.
Natural freeze tolerance in ectothermic vertebrates
.
Annual Review of Physiology
54
:
619
637
.
Storey,
K. B.,
and
J. M.
Storey.
1996
.
Natural freezing survival in animals
.
Annual Review of Ecology and Systematics
27
:
365
386
.
Tattersall,
G. J.,
and
R. G.
Boutilier.
1997
.
Balancing hypoxia and hypothermia in cold-submerged frogs
.
Journal of Experimental Biology
200
:
1031
1038
.
Todd,
J.,
J.
Amiel,
and
R.
Wassersug.
2009
.
Factors influencing the emergence of a northern population of eastern ribbon snakes (Thamnophis sauritus) from artificial hibernacula
.
Canadian Journal of Zoology
87
:
1221
1226
.
Ultsch,
G. R.
1989
.
Ecology and physiology of hibernation and overwintering among freshwater fishes, turtles, and snakes
.
Biological Reviews
64
:
435
516
.
Viitanen,
P.
1967
.
Hibernation and seasonal movements of the viper, Vipera berus berus (L.) in southern Finland
.
Pp
.
472
546
in
Annales Zoologici Fennici 4
.
Finnish Zoological and Botanical Publishing Board
,
Finland
.
Weatherhead,
P. J.,
and
K. A.
Prior.
1992
.
Preliminary observations of habitat use and movements of the eastern massasauga rattlesnake (Sistrurus c. catenatus)
.
Journal of Herpetology
447
452
.
Yagi,
A. R.,
and
R.
Tervo.
2005
.
[data sensitive] Massasauga population interim report
.
National Massasauga Recovery Team and the Ontario Ministry of Natural Resources
,
Vineland Station, Ontario, Canada
.
Yagi,
K. T.,
and
J. D.
Litzgus.
2012
.
The effects of flooding on the spatial ecology of spotted turtles (Clemmys guttata) in a partially mined peatland
.
Copeia
2
:
179
190
.
Zappalorti,
R. T.,
J.
Burger,
W. D.
Burkett,
D. W.
Schneider,
M. P.
McCort,
and
D. M.
Golden.
2014
.
Fidelity of northern pine snakes (Pituophis m. melanoleucus) to natural and artificial hibernation sites in the New Jersey pine barrens
.
Journal of Toxicology and Environmental Health Part A
77
:
1285
1291
.
Zuur,
A. F.,
E. N.
Ieno,
N. J.
Walker,
A. A.
Saveliev,
and
G. M.
Smith.
2009
.
Zero-truncated and zero-inflated models for count data
.
Pp
.
261
293
in M. Gail, K. Krickeberg, J. M. Samet, A. Tsiatis, and W. Wong (eds.),
Mixed Effects Models and Extensions in Ecology with R
.
Springer
,
USA
.

Appendix 1

Detailed Methods for Measuring the Life Zone (LZ).—Groundwater wells paired with frost tubes were installed, forming a grid pattern, across both Mined and Not-Mined hibernation study areas. Thirty-meter and 50-m spacings were used for Mined and Not-Mined well fields, respectively. Weekly measurements averaged 12 wells for Mined (0.8 ha) and 24 wells in the Not-Mined (3.75 ha) areas. We placed wells and frost tubes approximately 1 m away from any surface hole, to avoid impacting existing snake burrows. We made groundwater wells from 1.83-m length of black 5-cm diameter ABS pipe and drilled several 0.5-cm-diameter holes perpendicular through the sides of the bottom one-third portion, and then placed a cap on top (Fig. 3). To install a well, we bored a vertical hole into the ground surface up to 1.2-m depth using a hand-operated 5-cm-diameter soil auger. We placed the well into the hole and tamped down the soil around the pipe at the ground surface interface. We also measured the depth to clay (i.e., organic soil depth [cm]) at installation. Soil-depth measurement was limited to the maximum length of the auger ∼1.2 m. Fibric soil and subterranean spaces were also noted when they occurred.

We constructed a frost tube in two parts, following a design developed by the U.S. Forestry Service (Patric and Fridley, 1969). For the first part, we made an outer casing from 1.3-cm-diameter PVC conduit, and cut it into 60–100-cm lengths depending on soil depth. We covered one end with 3M™ duct tape and pushed the casing vertically into the ground, flush to the surface. We located the casing approximately 1 m from the well. For tighter soils, we recommend the use of metal rebar to make the hole. We constructed the inner second part from 0.95-cm outer-diameter clear, flexible tubing sealed at the bottom with a cured waterproof sealant and a 2-cm piece of wooden doweling to form a plug. The flexible tubing fits tightly inside the outer casing. We delineated the corresponding ground surface by marking the clear tube with an indelible felt-tipped pen. We filled the clear tube with red-food-dye colored water up to this mark. We used waterproof tape to secure the top of the frost tube to prevent water leakage and attached a guidewire through the tape to a nearby post or shrub to help locate the frost tube under the snow. We measured frost depth (FD [cm]) as ice formation within the red water-filled tube which changes color from red to clear, giving a distinct visual demarcation of ice. We measured GWL (cm) in the well from the top of the pipe using an electronic measuring tape (Heron Instruments Little Dipper™ 22 m), which beeps when it encounters water. LZ size (cm) equals the space between FD and GWL (Fig. 3). Manual data collection occurred weekly between 1 November and 31 April, from 2003 to 2014.

We also measured temperature (TGW [°C]) and dissolved oxygen (GWDO [mg/L]) within the top 10-cm groundwater layer inside the well, using a DO meter (YSI 55 or 550 models) calibrated to zero salinity, and surface elevation (183 m above sea level). We calculated the LZ groundwater thermal buffering function (TBUFF [°C]) for each study area as the difference between TGW and daily mean air temperature (TAIR; Environment and Climate Change Canada, 2016). We also measured total snow height (cm) including any ice accumulations, using a measuring tape. The flooded zone was mapped from field observations with the use of a handheld Garmin GPS (model eTrex 20x) and digitized into orthogonal corrected aerial imagery with ©ArcGIS mapping software. We recorded surface flooding and ice formation at a well site, and where this occurred, we scored LZ equal to zero centimeters.

Appendix 2

A summary of Massasauga population data collected during 11 active seasons following each winter's hibernation study (2004–2014) located in southern Ontario, Canada. Encounters (Nindiv), mean age (with standard error), detection rate (C+R / Nindiv), and flood state from 2004 to 2014 are included. Nindiv is the number of new captures, recaptures, and known-undetected snakes, calculated by backcasting age in the population time series. Flooding events started in the winter of 2006–2007 and continued through to winter of 2010–11.

A summary of Massasauga population data collected during 11 active seasons following each winter's hibernation study (2004–2014) located in southern Ontario, Canada. Encounters (Nindiv), mean age (with standard error), detection rate (C+R / Nindiv), and flood state from 2004 to 2014 are included. Nindiv is the number of new captures, recaptures, and known-undetected snakes, calculated by backcasting age in the population time series. Flooding events started in the winter of 2006–2007 and continued through to winter of 2010–11.
A summary of Massasauga population data collected during 11 active seasons following each winter's hibernation study (2004–2014) located in southern Ontario, Canada. Encounters (Nindiv), mean age (with standard error), detection rate (C+R / Nindiv), and flood state from 2004 to 2014 are included. Nindiv is the number of new captures, recaptures, and known-undetected snakes, calculated by backcasting age in the population time series. Flooding events started in the winter of 2006–2007 and continued through to winter of 2010–11.