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.
Materials and Methods
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.
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).
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.
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.
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.
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.
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).
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).
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).
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).
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.
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).
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.