Demographic Analysis of Imperiled Eastern Massasaugas (Sistrurus catenatus catenatus)

Precise estimates of demographic rates (e.g., survival and fecundity) are critical for the application of conservation tools like population viability analysis (PVA; Coulson et al., 2001). PVAs created with poor parameter estimates often result in inaccurate predictions of extinction rates and future population size (Beissinger and Westphal, 1998; Coulson et al., 2001; Ellner et al., 2002). Stochastic population models also require estimates of the process variance (e.g., temporal and spatial variance) of demographic rates. For example, a population model incorporating environmental stochasticity requires an estimate of the temporal process variance of survival as well as other demographic rates. Total variance is made up of the process variance(s) and the sampling variance. Because of the inclusion of sampling variance, using total variance in a population model causes an inflation of the probability of population decline and extinction (Gould and Nichols, 1998). In contrast, process variance accounts for the real variation in a parameter and is better suited for use in a PVA than the total variance (Gould and Nichols, 1998). Although population models for Eastern Massasaugas (Sistrurus catenatus catenatus) exist, models incorporating environmental stochasticity used the total variance of survival rather than the process variance (Middleton and Chu, 2004; Miller et al., 2005; Bissel, 2006; Bailey et al., 2011).

Other demographic characteristics like population size and growth rate also are useful parameters to estimate. Population size estimates are useful for determining trends in population size over the years of the study and whether the population is at a critically small size. Population growth rates also can be compared with projected growth rates from population models as a method to check the accuracy of the population model or to determine future trends in population growth (Miller et al., 2011).

Demographic estimates are needed for the conservation and management of imperiled species such as Eastern Massasaugas (Crother et al., 2012). Eastern Massasaugas are small, thick-bodied rattlesnakes found in the upper Midwest and Great Lakes regions of North America (Szymanski, 1998). Despite a relatively large range, populations are highly fragmented, and many populations appear to be in decline (Szymanski, 1998; Rouse and Wilson, 2002). Eastern Massasaugas are afforded some level of protection in every state or province in which the species resides (Szymanski, 1998). They are a candidate species for listing under the U.S. Endangered Species Act of 1973 (U.S. Fish and Wildlife Service, 2010) and are listed as threatened under Canada's Federal Species at Risk Act (Environment Canada, 2012). Because of this status, there is a demand for the application of tools like PVA and, therefore, estimates of demographic characteristics such as survival to inform management decisions.

Here, we use data from a long-term mark–recapture study of a northern population of Eastern Massasauga rattlesnakes to estimate a number of demographic parameters, including adult survival, sex differences in adult survival, temporal process variance of survival, population growth rate, and population size.

Data were collected from an Eastern Massasauga rattlesnake population on Beausoleil Island in southeastern Georgian Bay, Canada (Fig. 1). The island is ∼1,100 ha and has a moderate climate with reduced temperature extremes, delayed warming in spring, and milder fall temperatures compared to areas further inland. The southern two-thirds of the island consist of drumlins, till, and sand cover, whereas the northern portion consists of till and sand patches, ponds, streams, wetlands, and two small lakes. Vegetation also differs from north to south with mixed and coniferous forest in the north and deciduous stands in the south. Beausoleil Island has a long history of human occupation that includes multiple Native American tribes, the French, and the English. Parts of the island were farmed in the past, and some logging occurred; however, the island was used for recreation since the 1920s, and it was established as a national park in 1929. Currently, there are roads, trails, staff residences, and campgrounds on the island.

National park staff collected mark–recapture data from 1978–2008. The staff used three different marking methods over the years of the study: 1) from 1978–1985, plastic colored disks attached to rattle segments with monofilament; 2) from 1985–1992, subcaudal scale branding; and 3) from 1992–2008, passive integrated transponder (PIT) tags. Staff measured snakes for snout–vent length (SVL) and determined sex via probing the base of the tail at capture. Because these changing marking methods have different levels of efficacy, and to minimize the number of parameters in the models, we restricted the data analyzed to adult snakes in the 1992–2008 period. We defined adult snakes as animals with SVL > 350 mm based on a natural break in the size-frequency distribution observed in this population.

To estimate survival and recapture rates, we fit Cormack–Jolly–Seber (CJS) models to the data with the package “RMark” (Laake, 2013) for R 3.0.1 (R Core Team, 2014) and Program MARK 7.1 (White and Burnham, 1999). We fit the models using the log link function and second partial derivative method for variance estimation. We also estimated the process variance of the annual survival rate via variance components analysis of the most parsimonious model within MARK (Burnham et al., 1987; White et al., 2001).

The global model for this data set was the fully parameterized three-way interaction of sex, transience effect, and time for survival and the sex-by-time interaction for recapture rate. In terms of CJS models, transience is defined as individuals that are captured once but not recaptured throughout the rest of the study (Choquet et al., 2005). To account for these individuals, we included a transience effect (indicated by “tran” in the model parameters) in the candidate model set. The transience effect replaced typical time dependence in survival with survival over two different time periods: survival over the first year since marking and survival over subsequent years. We compared models using Akaike's Information Criterion corrected for small sample sizes (AIC_c). Our candidate set included all models nested within the global model (Table 1). We tested the global model for goodness-of-fit including testing for the transience effect (TEST3.SR and TEST3.Sm), testing for a trap dependence effect (TEST2.CT and TEST2.CL), and overall goodness-of-fit (sum of all tests) with U-CARE (Choquet et al., 2009). We used the median c-hat test within MARK to estimate overdispersion in the global model.

We generated population estimates using the Jolly–Seber model via the program JOLLY (http://www.mbr-pwr.usgs.gov/software.html; Pollack et al., 1990) with 95% confidence intervals calculated using Manly's method (Manly, 1984). We also calculated average population size over the study years via the Horvitz–Thompson estimator (Horvitz and Thompson, 1952). We estimated the population growth rate with the Pradel survival and population growth rate model within MARK (Pradel, 1996). The same model fitting procedure used for the CJS models above was used for the Pradel models.

The mark–recapture data set consisted of 262 marked snakes (160 males and 102 females) captured 431 times from 1992–2008. The sum of all tests with pooled sexes indicated no sex effect (χ² = 45.43, df = 55, P = 0.82); thus, all further goodness-of-fit tests were performed on the pooled sex data set. The test for transience was significant (χ² = 30.31, df = 15, P = 0.01); the test for trap-dependence was non-significant (χ² = 7.23, df = 14, P = 0.92). Seventeen models were fit to the data (Table 1). The most parsimonious (model weight = 0.79) contained the parameters “tran” for survival rate and “sex” for recapture rate. The second most parsimonious model (model weight of 0.20) contained the parameters “sex” and “tran” in interaction for survival rate and “sex” for recapture rate. We used model averaging to calculate survival and recapture estimates (Table 2). The survival rate estimate for the first year since marking was lower than the survival rate estimate in subsequent years for both males and females, but there were only small differences in magnitude between males and females for either survival rate estimate (Table 2). Recapture rates were consistent across years, but differed between the sexes (Table 2).

The adult survival model “Phi(sex × tran × time) p(time)” (sex, age, and time dependent survival and time dependent recapture rates) was used to calculate the process variance of the survival rate within MARK because time-dependence is necessary for the estimation process (White et al., 2001). The process variance for males was 0.006 and accounted for ∼11% of the total variance (0.054). The process variance for females was unestimable, perhaps attributable to insufficient data.

Population size estimates from Program JOLLY ranged from 35–77 (mean = 55) over the study years (Fig. 2). The Horvitz–Thompson estimate for average population size over the study years was 54 ± 4.03 SE. Two Pradel survival and population growth models based on the most parsimonious CJS model (Phi[tran] p[sex]) were fit within MARK: 1) with time-dependent population growth rates and 2) with time-independent population growth rate. The model with the time-independent population growth rate was clearly the most parsimonious (Table 3). The growth rate estimate was 1.02 (95% CI = 0.99–1.04), indicating a stable to slightly increasing population.

This study provides an unusually complete analysis of survival and population size of an Eastern Massasauga population near the northern limits of the species' range. Importantly, we demonstrate a significant effect of transience on adult survival. The models including the transience effect contained two separate survival rates that varied greatly from each other. The first survival rate (0.44 for males; 0.47 for females) is the apparent survival of snakes during their first year after marking; the second survival estimate (0.74 for males; 0.73 for females) is the survival during subsequent years (i.e., of resident snakes). Apparent survival is the product of the actual survival probability and probability of transience. The transience effect in the model effectively partitioned the effect of transient snakes in the apparent survival rate of the first year since marking. Because we are interested only in the probability of survival over the bulk of the study, the apparent survival rate of the first year since marking can be considered a nuisance parameter. We consider the second survival estimate to reflect more accurately the true survival rate. One biological interpretation of this transience effect is that two subclasses of snakes are present; a “resident” subclass that was sampled consistently throughout the study and a “transient” subclass consisting of individuals that were sampled only rarely (e.g., as they moved onto and off of Beausoleil Island or as they moved in and out of more consistently sampled parts of the island). Whether such transience characterizes Eastern Massasaugas or other snake species more generally warrants investigation. Although this study considers the Eastern Massasaugas residing on Beausoleil Island as a single population, there is movement to and from the island (K. Prior, pers. comm.). The population on Beausoleil Island is part of a larger, roughly defined metapopulation found on the eastern side of Georgian Bay (Szymanski, 1998; Parent and Weatherhead, 2000).

Previously, Middleton and Chu (2004) attempted to use this same data set to estimate a survival rate for a PVA; however, the spreadsheet containing the data was corrupted (possibly attributable to a sorting error) prior to their analysis (Middleton and Chu, 2004). They were able to extract some demographic data and perform a time series analysis on the data set, but they did not estimate survivorship (Middleton and Chu, 2004). Instead, they based their survivorship estimate (0.69) on information from multiple sources, including components of the data set analyzed in the current study (Middleton and Chu, 2004). Brennan and Tischendorf (2004) and Miller et al. (2005) based their model parameters on Middleton and Chu's work and used an adult survival rate of 0.69 as well. The annual survival rate estimated in the current study was somewhat higher than the rates used in those Eastern Massasauga PVAs, although the difference between 0.74 and 0.69 may not be significant. Other adult survival estimates for Eastern Massasauga populations range from 0.33–0.95 (table 1 in Jones et al., 2012) and encompass our estimates for both males and females. Adult survival estimates generated using known-fate models applied to radiotelemetry data show a geographic trend of increasing survival from the southwest to the northeast (Jones et al., 2012). Our average adult annual survival estimate was at the higher end of the range of survival estimates for other Eastern Massasauga populations in Ontario (0.62–0.79) and was higher than the average survival across populations of 0.67 (Jones et al., 2012).

Calculation of the process variance requires long-term data sets. Few other estimates of the process variance of survival have been published for snakes or other reptiles (Whiting et al., 2008; Miller et al., 2011; Stanford, 2012). Process variance of adult and subadult survival for four populations of Nerodia sipedon insularum ranged 0.0005–0.125 (Stanford, 2012). Process variance for adult survival in Nerodia harteri paucimaculata was 0.0073 over 9 yr (Whiting et al., 2008). Process variance of adult survival in five populations of Thamnophis elegans range 0.00003–0.00353 over 8–13 yr (Miller et al., 2011). We were unable to estimate process variance for females; thus, we focus on the male estimate here. The estimate of 0.006 for the male Georgian Bay population falls within the process variance range for N. s. insularum, indicating similar variation in survival over time but was much higher than estimates for T. elegans or N. h. paucimaculata populations, indicating much more variation in survival over time. The unrealistic estimate of process variance for females may be attributable to the smaller sample size of females.

Previous Eastern Massasauga PVAs used the total variance of population parameters (table 1 in Jones et al., 2012). Total variance is made up of the process (= temporal) variance and the sampling (= error) variance, and its use in a PVA will result in an inflated estimate of extinction risk (Gould and Nichols, 1998). As a demonstration of this inflation, we created two population models in VORTEX 10.0.7.9 (Lacy and Pollak, 2014) based on Eastern Massasauga demographic values (Parent and Weatherhead, 2000; Dreslik et al., 2011), including the adult survival and variance estimates from this study. Parameters were adjusted to give the models a deterministic population growth rate of 1.00. The models were identically parameterized except we used the total variance of adult survival as the measure of environmental stochasticity in the first and the process variance of adult male survival in the second (variance was set to 0 for all other survival parameters). Over 1,000 simulations, the probability of extinction over a 50-yr period was 20.8% for the model using total variance, but the probability of extinction was only 3.0% for the model using process variance.

The mean population size estimates of 55 and 54 are quite small, and accordingly, this population is at risk of extinction (Seigel and Sheil, 1999; Mauger and Wilson, 1999). Although our estimation of population size did not use a maximum likelihood methodology, the relative stability of our population size estimates over the years is supported by the population growth rate estimate from the Pradel model. Genetic data suggest Eastern Massasauga populations are frequently small (100 to a few thousand individuals) and relatively isolated (Chiucchi and Gibbs, 2010). In a small population such as this, a 1–2% increase in annual adult mortality can cause large decreases in population size and ultimately population extinction (Seigel and Sheil, 1998). Given the Eastern Massasaugas' protected status in Canada, however, the population may remain stable or even increase in size if protection continues, because this population resides on government land and our long-term population growth rate estimate is slightly greater than 1.00.

Reconstruction and analysis of this 30-yr Eastern Massasauga data set produced multiple demographic values, including the first estimate of temporal process variance in survival for this species and one of the few for any snake species. These demographic values are useful not only for assessing the health of this population of Eastern Massasaugas but also for comparison with known values of other populations and in future analyses such as population projection models and PVAs. As for the Beausoleil Island population itself, retrospective analysis showed the population was relatively stable for over a decade. If the population remains undisturbed, this stability may continue. Further monitoring and conservation-related management will help ensure this Eastern Massasauga population's persistence.

We thank M. Desjardins, K. Prior, and A. Promaine for their assistance in obtaining the dataset analyzed here and reviewing the manuscript. We also thank the park staff who collected much of the mark-recapture data, and C. Jaeger for assistance with the map.