Eastern Indigo Snakes (EIS, Drymarchon couperi) and Eastern Diamondbacked Rattlesnakes (EDB, Crotalus adamanteus) are species of conservation concern, in large part attributable to anthropogenic landscape changes within the southeastern Coastal Plain of North America. Both species use Gopher Tortoise (Gopherus polyphemus) burrows on xeric sandhills for winter retreat sites. Protected lands play an important role in the conservation of threatened species by offering the potential to conserve potentially limiting resources such as sandhills. We surveyed 40 randomly selected xeric sandhills containing Gopher Tortoise burrows on protected lands throughout the Lower Altamaha River Watershed in southern Georgia using visual encounter surveys over three winters (November through March). We used single-season occupancy models to relate detection and occupancy rates to survey- and site-specific covariates collected at both the sandhill- and landscape-scale. Eastern Indigo Snake occupancy was positively related to the number of Gopher Tortoise burrows and the amount of surrounding sandhill habitat. In contrast, EDB occupancy was not associated with any of the covariates we considered, perhaps because EDB/EIS use a greater diversity of winter retreat sites. Detection of EIS was higher than EDB (0.40 vs. 0.22) and most influenced by air temperature, whereas EDB detection was most influenced by survey date. Our study provides previously lacking population-level detection rates and habitat associations for EIS and corroborates the previously noted importance of Gopher Tortoise burrows as overwintering retreat sites. Our study also illustrates the potential shortcomings of monitoring multiple species using survey methodologies designed for a single species.
Anthropogenic alterations to natural landscapes represent one of the most severe threats to biodiversity (Wilson, 1999; Noss et al., 2006). Because many imperiled species cannot persist in human-modified habitats, protecting areas of natural habitat is essential for their conservation. As a result, despite the controversies and complexities of establishing wildlife habitat reserves, protected lands remain a vital component of conservation efforts (Bruner et al., 2001; Ervin 2003; Watson et al., 2014). Beyond simply protecting natural habitats, however, management actions (including assisted colonization and increased law enforcement) that incorporate greater knowledge of species' life histories may be required to prevent further population declines, even on protected lands (Carroll et al., 2004; Liu et al., 2001; McCoy et al., 2006; Turner et al., 2006; Williams et al., 2012). Populations on conservation lands also may suffer the negative effects of population isolation and demographic stochasticity if population connectivity is not maintained (Cushman et al., 2012; Trainor et al., 2013).
The longleaf pine (Pinus palustris) ecosystem (LLPE) formerly dominated much of the Coastal Plain of the southeastern United States (Landers et al., 1995; Van Lear et al., 2005), but urban development, agriculture, and modern silviculture practices have eliminated ∼97% of the original habitat (Van Lear et al., 2005). ). The LLPE formerly supported regular low-intensity fires during the spring and summer, which maintained an open, savannah-like structure and a diverse ground cover of grasses and forbs (Waldrop et al., 1992; Platt, 1999; Van Lear et al., 2005). Widespread fire suppression has resulted in substantial structural changes to much remaining habitat by increasing hardwood densities (e.g., oak, Quercus spp.) that, in turn, increased canopy cover and reduced native ground cover species (Engstrom et al., 1984; Gilliam and Platt, 1999). As a result of these changes, many LLPE-associated taxa have declined (Guyer and Bailey, 1993; Trani-Griep, 2002; Means, 2006). Eastern Indigo Snakes (EIS, Drymarchon couperi) and Eastern Diamondbacked Rattlesnakes (EDB, Crotalus adamanteus, Crother, 2012) are closely associated with the LLPE throughout much of their ranges (Means, 2006; Waldron et al., 2008; Enge et al., 2013), although both species also occur in non–longleaf pine-dominated habitats, particularly in peninsular Florida (e.g., scrub and maritime forests; Moler, 1992; Jensen et al., 2008). Both species also have undergone declines throughout their ranges, particularly the EIS that is listed as federally threatened (United States Fish and Wildlife Service, 1978; Martin and Means, 2000; Timmerman and Martin, 2003).
Both species have several life-history traits, which have likely exacerbated the effects of anthropogenic habitat and land-use changes. First, both species have large spatial requirements, particularly the EIS whose annual home ranges can exceed 1,500 ha (Waldron et al., 2006; Breininger et al., 2011; Hoss et al., 2010; Hyslop et al., 2014). Second, both species use a diversity of habitats and may use habitats according to season (Hoss et al., 2010; Hyslop et al., 2014). Third, within sympatric parts of their northern distributions (i.e., southern Georgia and northern Florida), both species use Gopher Tortoise (Gopherus polyphemus) burrows as winter (i.e., cool-season, November through March) retreat sites and, hence, are closely associated with xeric sandhills where tortoise burrows may be found. The EDB, whose range extends northward beyond the range of the Gopher Tortoise, also will use alternate winter retreat sites (e.g., stump holes, root channels, Nine-Banded Armadillo [Dasypus novemcinctus] burrows; Martin and Means, 2000; Timmerman and Martin, 2003), whereas the EIS uses tortoise burrows almost exclusively within the northern part of its distribution (Diemer and Speake, 1983; Stevenson et al., 2003; Hyslop et al., 2009). Finally, both species show high interannual fidelity to overwintering sites (Stevenson et al., 2003, 2009; Waldron et al., 2013a) that may prove maladaptive if habitat changes render their overwintering habitat unsuitable. Given these traits, protected lands within the LLPE may greatly benefit EIS and EDB through the preservation of potential overwintering habitat.
The Lower Altamaha River Watershed (LARW) contains an important portion of EIS and EDB distributions within the state of Georgia (Martin and Means, 2000; Enge et al., 2013). Although protected lands are widespread within this watershed, many have been highly altered due to fire suppression or commercial forestry practices. Determining patterns of winter sandhill occupancy by EIS and EDB within these lands is important for determining their conservation value to those species. Our overall goal was to determine EIS and EDB occupancy of potential overwintering habitat on conservation lands within the LARW and identify factors influencing occupancy rates. Our specific objectives were to 1) estimate and compare EIS and EDB occupancy of potential overwintering habitat within protected lands across our study area; 2) evaluate hypothesized relationships between detection rates and survey-level covariates; and 3) evaluate hypothesized relationships between winter occupancy of xeric sandhills and sandhill- and landscape-scale covariates.
Materials and Methods
Study Area.—Our study area was LARW within the current range of the EIS (Enge et al., 2013; Fig. 1). We excluded the Oconee River watershed because contemporary EIS records are lacking for this drainage (Enge et al., 2013). Natural upland communities included xeric sandhills, pine flatwoods, mixed oak-pine forests, and slope forests (i.e., oak, beech [Fagus spp.], and magnolia [Magnolia spp.] communities). Historically, sandhills and pine flatwoods were dominated by longleaf pine, but many were since converted into commercial pine forests of slash (Pinus elliottii), loblolly (Pinus taeda), or sand pine (Pinus clausa). Widespread and characteristic wetland habitats included bottomland hardwood forests, bay swamps, blackwater river and floodplain swamps, creek swamps, and seasonal depressional ponds forested with cypress (Taxodium spp.) and/or gum (Nyssa spp.).
Fig. 1.
Map of study area and Eastern Indigo Snake (EIS) and Eastern Diamondbacked Rattlesnake (EDB) occupancy monitoring sites by subdrainage within the Lower Altamaha River Watershed (LARW) in southeastern Georgia. The contemporary distribution of the EIS and EDB are shown in the insert in dark and light grey, respectively, following Enge et al. (2013) and Jensen et al. (2008). Boundary lines within the LARW denote watershed boundaries used to stratify our sampling of survey sites.
Site Selection.—We selected several state- and privately owned properties within the LARW that were recognized for their conservation value to the LLPE or its associated biodiversity and that we felt reasonably certain we could access indefinitely as part of our ongoing monitoring program. We initially conceived this to be an EIS monitoring study; because EIS almost exclusively use tortoise burrows for winter retreat sites (e.g., Hyslop et al., 2009) and the value of protected lands for overwintering habitat is poorly understood, we identified discrete xeric sandhill patches as spatial sampling units (“sites”). All potential sites were identified as unique sandhill polygons from the Georgia Department of Natural Resources (GA DNR) Nongame Conservation Section's sandhills GIS layer (Elliott, 2009). We refined our pool of potential sites to those supporting multiple (>1) tortoise burrows that we identified using previous observations or a combination of visual assessments of aerial imagery and ground truthing. After identifying all potential sites within accessible properties, we randomly selected 34 sites that were surveyed in 2010–2013. We added 6 additional randomly selected sites in the second year such that 40 sites were surveyed in 2011–2013. To ensure the sites were distributed throughout our study area, we stratified our sampling by watershed and property type (privately vs. government or nonprofit owned) with a greater emphasis (58% of sites) on government or nonprofit owned properties because of their perceived importance for EIS conservation. Mean (± SD) site size was 35.16 ± 42.96 ha (range = 3.17–215.13 ha), and the mean distance between a site's centroid and the centroid of its nearest neighboring site was 2.79 ± 3.76 km (range = 0.42–19.26 km).
Survey Procedures.—We visited each site four times between 1 November and 31 March (i.e., in winter) in each year from 2010 to 2013, yielding three winter sampling seasons. At sites that were too large to survey within a single day, we randomly selected one corner of the sandhill as a starting point, surveyed as much of the sandhill as possible, and considered the area surveyed as the extent of the site. We conducted visual encounter surveys with one to three observers inspecting tortoise burrow entrances for snakes or their shed skins (Stevenson et al., 2003, 2009). On the first site survey, observers marked as many tortoise burrows (active, inactive, and abandoned; Auffenberg and Franz, 1982) as possible using GPS units, and used these marked points to guide subsequent survey efforts. Because tortoise burrows were typically not distributed throughout the entire sandhill, the actual area surveyed often was smaller than the size of the site. The vast majority of surveys (95%) were conducted by a single observer. If we observed a fresh and distinct snake track at a burrow, we scoped the burrow using a burrow camera system (CCD Hi Resolution Black and White Camera and Black and White Active Matrix Backlit Widescreen LCD Monitor; Sony, Inc., Tokyo, Japan) at the end of the survey to confirm the presence of an EIS or EDB. Surveys were conducted between 0900 and 1700 h, and we did not conduct surveys if the forecasted air temperature high for the day was <10°C.
Occupancy Modeling.—We modeled our data using single-season occupancy models (MacKenzie et al., 2002). Because EIS and EDB are long-lived species (Bowler 1977; Stevenson et al., 2003; Waldron et al., 2013a) that exhibit high fidelity to overwintering sites (Stevenson et al., 2003; Waldron et al., 2013a), turnover rates (colonization and extinction) were of less interest than the determinants of occupancy; and as such, we conducted a single-season, rather than multiseason, analysis. Year-to-year variation in site occupancy can easily be accommodated by “stacking” the data such that each site-year combination is represented as a unique site (4 visits at 114 sites) and enforcing an additive “year” effect on occupancy (e.g., Miller et al., 2013).
Hypotheses for Occupancy Rate.—We hypothesized that EIS and EDB winter occupancy of xeric sandhills would be associated with multiple factors operating at both the scale of the sandhill and the surrounding landscape (Table 1). Tortoise burrow abundance was calculated as the mean number of tortoise burrows across all surveys of each site. We estimated the amount of potential overwintering habitat at each site by manually digitizing GIS polygons containing the GPS-marked tortoise burrows at each site and calculating the area of those polygons. Tortoise burrow abundance and area of potential habitat were positively correlated (rs = 0.59) and, therefore, never were included in the same model. We also calculated an index of tortoise burrow density by dividing the number of tortoise burrows by the area of potential habitat.
Table 1.
Site- and landscape-scale covariates and their hypothesized relationships between Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) winter occupancy of xeric sandhills supporting Gopher Tortoise (Gopherus polyphemus) burrows in the Lower Altamaha Watershed in southern Georgia. We provide the hypothesized relationship (±) and relevant references. See text for details of covariate measurements.
We hypothesized that vegetation structure would influence occupancy because of its influence on Gopher Tortoise habitat suitability (Aresco and Guyer, 1999; Boglioli et al., 2000; Bauder et al., 2014; Kowal et al., 2014) and potential for creating suitable microhabitats for thermoregulation (Rubio and Carrascal, 1994; Blouin-Demers and Weatherhead, 2002). We measured vegetation variables at randomly selected points at each site proportional to the size of the site (14–30 points per site). We used the point-center-quarter method (Cottam and Curtis, 1956; Beasom and Haucke, 1975) to record the distance from each sampling point to the nearest tree (defined as ≥2.5 cm DBH) within four quadrants formed by the cardinal directions (i.e., NE, NW, SW, and SE). For each nearest tree, we measured the tree's diameter at breast height (DBH) and classified tree species as pine, oak, and other hardwood. We then calculated relative density, relative basal area, and relative frequency for each species group and summed these values across sampling points to generate importance values for each species group at each site (Cottam and Curtis, 1956). We dropped the importance value of other hardwoods because of insufficient data. The importance values for oaks and pines were highly correlated (rs = −0.98); thus, we retained the importance value of pine. We recorded canopy cover using a spherical densiometer at points 5 m from the sampling point in each of the cardinal directions and calculated the mean and SD of canopy cover for each site across all sampling points. We estimated shrub cover within a 10 m radius around each sampling point as one of five categories (0%, 1–25%, 26–50%, 51–75%, 76–100%) and converted these categories into ordinal variables that were used to calculate the mean and SD of shrub cover for each site.
We recorded sandhill condition using the sandhills GIS layer developed by GA DNR (Elliott, 2009). Condition was a measure of ecological integrity and was subjectively determined based on soil disturbance, vegetation density, and presence of nearby development and other anthropogenic encroachment visible from aerial imagery. Sandhill condition was recorded as excellent (N = 1), good (N = 9), fair (N = 24), and poor (N = 6); we combined the excellent and good categories. We hypothesized that EIS presence would negatively affect the presence of EDB, because EIS will predate EDB (Stevenson et al., 2010; Steen et al., 2014). Because our sample sizes were too small to consider a multispecies model (MacKenzie et al., 2004b), we considered a binary covariate denoting whether EIS was detected at a site at any point during our 3-yr study.
We measured landscape variables using buffers of varying radii centered on each site to test for multiscale occupancy-covariate associations (Johnson et al., 2004; Table 1). We used published estimates of EIS and EDB year-round home-range sizes (Waldron et al., 2006; Hoss et al., 2010; Hyslop et al., 2014) to select buffers whose sizes approximated the mean and maximum home-range sizes of each species. These sizes ranged from 359–1,530 ha for EIS and 29–62 ha for EDB; however, species may respond to landscape features at scales beyond their home range (Kie et al., 2002). For example, Steen et al. (2012a) found that EDB occupancy was related to landscape composition within 315-ha buffers. Furthermore, male EIS may make linear movements from overwintering sites to summer foraging habitats 1.5–7.5 km in length (Hyslop et al., 2014). Therefore, we used the following ranges of buffer radii for each species: 0.25, 0.50, and 1.00 km for EDB and 1.00, 2.00, and 5.00 km for EIS, which resulted in 19.63-, 78.54-, 314.16-, 1,256.64-, and 7,853.98-ha buffers, respectively.
We used the Georgia Land Use Trends (GLUT) 2008 land cover data (Natural Resources Spatial Analysis Laboratory, University of Georgia, Athens) to measure the amount of clearcut/sparse forest, deciduous forest, evergreen forest, mixed forest, and wetland (forested and non-forested wetlands combined) cover types within each buffer. The amount of deciduous forest cover type was correlated with the amount of mixed forest cover type across all buffer sizes (rs = 0.53–0.87); therefore, we retained only mixed forest. To measure heterogeneity in landscape configuration, we calculated the number of habitat patches using wetland, clearcut, mixed forest, deciduous forest, and evergreen forest cover types (Couturier et al., 2014) and the edge density of those cover types within each buffer using the “SDMTools” package (v1.1-221, VanDerWal et al., 2014) in R v3.0.2 (R Foundation for Statistical Computing, Vienna, Austria). We also measured the amount of sandhill (regardless of condition) within each buffer using GA DNR's sandhill GIS layer (Elliott, 2009). We measured the amount of impervious surface within each buffer using the GLUT 2008 impervious surface raster layer. Impervious surfaces included unpaved roads which could have a negative impact on EIS and EDB through direct road mortality and increased human access to remote areas that could lead to increased persecution (Shepard et al., 2008; Robson and Blouin-Demers 2013; Breininger et al., 2012; but see Steen et al., 2007). Finally, we measured the amount of agricultural land cover using the GLUT land cover data. All continuous site covariates were z-score standardized prior to analysis with mean = 0 with standard deviation = 1 to facilitate model convergence (MacKenzie et al., 2006).
Hypotheses for Detection Probability.—To evaluate a priori hypotheses about variation in detectability, we measured specific survey-level covariates. Because males of both species actively search for females during the breeding season, we hypothesized that detection would be highest during the peak breeding season (November through January for EIS and August through October for EDB; Stevenson et al., 2009; Hoss et al., 2011; Waldron et al., 2013a). Winter surface activity also may be associated with higher air temperatures (Spence-Bailey et al., 2010; Couturier et al., 2013); therefore, we considered both a linear effect of air temperature and a quadratic effect of survey date that we considered a proxy for air temperature. We hypothesized that observer variability could affect detection rate (Lotz and Allen 2007; Alldredge et al., 2007); thus, we included a three-level categorical variable to denote which observer conducted the survey. We hypothesized that detection would increase with the amount of area surveyed (Chen et al., 2009); hence, we considered a linear effect of amount of potential habitat. Finally, we hypothesized that previous knowledge of species' detection at a site could bias survey efforts (Riddle et al., 2010), particularly because randomly assigning observers among sites was logistically infeasible. Therefore, we considered two binary covariates that denoted whether a species had been previously detected at a site (Riddle et al., 2010). One covariate denoted whether a species was detected previously at a site during a given winter (e.g., November 2010 through March 2011) and the other at any point during the study (MacKenzie et al., 2004a). For example, a detection history of 0010-0000-0101 (detected on the third survey in the first winter, not detected during the second winter, and detected on the second and fourth survey of the third winter) would have sampling covariates for that site of 0001-0000-0011 for detected during a given winter and 0001-1111-1111 for detected during the study (MacKenzie et al., 2004a). We z-score standardized all continuous sampling covariates.
Model Selection.—Prior to model fitting, we examined our final set of covariates for colinearity and did not include correlated (rs >0.60) covariates within the same model. We fit our models using the package “unmarked” (v0.10-4, Fiske and Chandler, 2011). We estimated the overdispersion parameter (c-hat) for each species using the global model, defined as the most parameter-rich model for both detection and occupancy that would converge, using the parametric bootstrap procedure of MacKenzie and Bailey (2004) in the “AICcmodavg” package (v2.0-3, Mazerolle, 2015). If c-hat was >1, we evaluated models using AIC adjusted for overdispersion (QAIC, Burnham and Anderson 2002). We first selected the best-supported detection model by using the most parameter-rich occupancy model and evaluating competing detection models. We retained detection covariates from models with ΔAIC/QAIC ≤ 2 for all subsequent analyses. We initially analyzed live snakes and shed skin separately, but our results were similar; therefore, we report the results of the pooled detections. To select the appropriate scale for each landscape-scale covariate, we conducted a stepwise procedure where, for each landscape-scale covariate, we compared that covariate measured across each set of scales for each species using AIC/QAIC and selected the best supported scale for inclusion in the final model set. We drew inference from detection and occupancy covariates in models with ΔAIC/QAIC ≤ 2 and considered a covariate significant if its 95% CI did not include zero. We calculated predicted values for plotting relationships between detection/occupancy and covariates using the modavgPred function (Mazerolle, 2015). We report derived estimates of detection and occupancy for each species as the model-averaged predicted value and 95% CI for detection and occupancy holding all sampling and site covariates constant at their mean values. Because our sites were clustered in space, we tested for residual spatial autocorrelation (rSAC) by plotting correlograms of Moran's I calculated on the residuals (Moore and Swihart, 2005). If we detected significant (P > 0.05) autocorrelation, we included an autocovariate term in the model to account for rSAC (Augustin et al., 1996; Moore and Swihart, 2005). We then reran our analyses including autocovariate terms in models where warranted. Additional details of our spatial autocorrelation analyses are provided in Appendix 1.
Results
Eastern Indigo Snakes.—We detected EIS at 50 of 114 (0.44) “site-years” over our 3-yr study. Live snakes and shed skins were detected during 8% and 6% of surveys, respectively. We did not detect a consistent, significant pattern of rSAC in our best ranked model (the “global” site- and landscape-scale model containing number of tortoise burrows; Appendix 2) but found significant rSAC within our other top-ranked models at multiple distance bins. Although adding autocovariates to these models dampened rSAC, we were unable to remove the significant rSAC within the 0–800 m bin (e.g., Appendix 2). This suggests some non-independence among closely spaced sites, consistent with our model overdispersion (c-hat = 1.61; Lebreton et al., 1992; MacKenzie and Bailey, 2004). Although this is unlikely to bias parameter estimates, it may inflate SE's and increase the risk of committing Type I errors (Griffith, 2003; Haining et al., 2009). We suggest these effects are minimized by adjusting our standard errors by the square root of c-hat, which then widens the 95% CI. Generally, our sites were ecologically independent (Appendix 1), because none of the distances between nearest neighbors, measured using site centroids, were less than the maximum diameter of reported EIS winter home-range sizes, and only 35% were less than the maximum diameter when distance between sites was measured in relation to site edges (Appendix 3).
Models containing air temperature and whether a detection was made previously within a winter had the best support (i.e., ΔQAIC < 2; Table 2). We retained air temperature in the detection model for all subsequent analyses although our inferences regarding occupancy were similar when we used our binary covariate denoting a previous detection within a winter to model detection. The model-averaged beta estimate and 95% CI for air temperature indicated a significant increase in EIS detection with increasing air temperature (β = 0.68, 0.05–1.31; Fig. 2A). The model-averaged beta estimate for previous detection within a winter suggested that EIS detection decreased once EIS was detected at a survey site, although this effect was not significant (β = −1.13, −2.29–0.02).
Table 2.
Model selection results, parameter estimates (β), and 95% CI for Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) detection during winter surveys of xeric sandhills supporting Gopher Tortoise (Gopherus polyphemus) burrows. Occupancy was modeled using the most parameter-rich occupancy model from our candidate set of occupancy models. Models are ranked according to Akaike's Information Criteria adjusted for overdispersion (QAIC). The overdispersion parameter (c-hat) was 1.83 for EIS and 1.00 for EDB. Deviance (Dev) is calculated as −2*quasi-log-likelihood, K represents the number of parameters in the model, and wi is the model weight. Multiple rows show the parameter estimates and CI for models with quadratic effects or multilevel categorical covariates.
Fig. 2.
Relationships between Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) probabilities of detection (p) and occupancy (psi) and covariates from the best supported (ΔQAIC ≤ 2) models with 95% CI that did not include zero: (A) air temperature (°C); (B) survey date (1 = 1 November); (C) number of Gopher Tortoise (Gopherus polyphemus) burrows; (D) amount of sandhill within the 1-km buffer; and (E) importance value of pine (Pinus spp.). Solid lines represent the model-averaged predicted values and the gray shaded band represents the 95% confidence interval.
Landscape covariates measured using 1 or 5 km buffers generally had the strongest support among spatial scales for EIS, but model uncertainty was relatively high (ΔQAIC ≤ 3.66; Appendix 4). Three of our 20 final candidate models had ΔQAIC ≤ 2, the “global” site- and landscape-scale combination model containing number of tortoise burrows, the model containing only number of tortoise burrows, and the model containing only the amount of sandhill within the 1-km buffer (Table 3). EIS occupancy was significantly associated with increasing numbers of tortoise burrows and amount of sandhill within 1-km (Table 4; Figs. 2C, 2D). The combination scale model also indicated that EIS occupancy was significantly associated with increasing pine importance (Table 4; Fig. 2E). The model containing only pine importance had low model support (wi = 0.02), however, suggesting the significant association with pine importance is conditional upon increasing numbers of tortoise burrows.
Table 3.
Model selection results, parameter estimates (β), and 95% CI for Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) occupancy during winter surveys of xeric sandhills supporting Gopher Tortoise (Gopherus polyphemus) burrows. Additive effects of year and, for EIS, a spatial autocovariate term were included in all models. Detection was modeled using the QAIC-best (AIC adjusted for overdispersion) covariate for detection. The overdispersion parameter (c-hat) was 1.61 for EIS and 1.00 for EDB. Deviance (Dev) is calculated as −2*quasi-log-likelihood, K is the number of parameters in the model, Δ is the ΔQAIC, and wi is model weight. See Table 1 for descriptions of model covariates.
Table 4.
Parameter estimates (β) and 95% CI for occupancy covariates for Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) occupancy during winter surveys of xeric sandhills supporting Gopher Tortoise (Gopherus polyphemus) burrows. Only results from models with ΔAIC/QAIC ≤ 2 are presented. Standard errors were adjusted by the square-root of the overdispersion parameter (c-hat = 1.61 and 1.00, for EIS and EDB, respectively). See Table 1 for a description of the covariates. AC is the autocovariate term used to model residual spatial autocorrelation and the distance refers to the neighborhood size. Inverse Euclidean distance was used to calculate AC.
Our derived estimate (i.e., model-averaged prediction with all covariates at their mean value) for EIS detection was 0.40 (95% CI = 0.27–0.55). Our derived estimates for EIS occupancy were 0.33 (95% CI = 0.14–0.60) for year 1, 0.19 (95% CI = 0.07–0.43) for year 2, and 0.27 (95% CI 0.11–0.53) for year 3.
Eastern Diamondbacked Rattlesnakes.—We detected EDB at 47 of 114 (0.41) “site-years” over our 3-yr study. We detected EDB shed skins during only one survey where live EDB were not detected. Overdispersion in the global model was low (c-hat = 1); therefore; we used AIC to evaluate models. We found relatively little rSAC in our EDB models and adding autocovariate terms did not improve model fit or reduce rSAC (Appendix 2); hence, we did not include autocovariate terms in our final analysis. Models containing a linear effect of survey date, a quadratic effect of survey date, and an observer effect on detection received the best support (max. ΔAIC ≤ 1.50; Table 2). Only the 95% CI for the linear effect of survey date did not include zero and indicated that detection decreased as the survey season progressed (Fig. 2B). We retained the linear effect of date, although our inferences regarding occupancy were robust to the detection model we used.
Four of our nine landscape scale covariates had the highest model support when measured within the 0.25-km buffer, although model uncertainty was relatively high across all landscape covariates for EDB (max. ΔAIC ≤ 2.97; Appendix 5). Six models had strong support (ΔAIC ≤ 2), although uncertainty was relatively high across these models (max. w = 0.17, Table 3). All of the 95% CI for the beta estimates in these models included zero, however, although the CI for evergreen forest and importance value of pine had a very small amount of overlap with zero, suggesting that EDB occupancy is negatively associated with the amount of evergreen forest and pine importance (Table 4).
Our derived estimate for EDB detection was 0.22 (95% CI = 0.14–0.34). Our derived estimates for EDB occupancy were 0.39 (95% CI = 0.15–0.69) for year 1, 0.36 (95% CI = 0.17–0.61) for year 2, and 0.59 (95% CI 0.29–0.84) for year 3.
Discussion
Our study indicates that both site- and landscape-scale covariates influenced EIS winter occupancy of xeric sandhills supporting tortoise burrows. In contrast, our covariates had little influence on EDB winter occupancy of these sandhills. We suspect this difference reflects the broader use of winter retreat sites and/or habitats by EDB compared to EIS. Although EIS appear to show near-exclusive use of Gopher Tortoise burrows on xeric sandhills for winter retreat sites (Stevenson et al., 2003, 2009; Hyslop et al., 2009, 2014), EDB will also use root holes, stump channels, and mammal burrows (Martin and Means, 2000; Timmerman and Martin, 2003). Furthermore, EDB are found northward beyond the range of the Gopher Tortoise within the southeastern Coastal Plain (Martin and Means, 2000). Our study was conceived and designed with an emphasis on monitoring winter EIS occupancy on protected lands and, therefore, our site selection and covariates focused on potential EIS overwintering habitat. As a result, our study sites may have captured some unknown, but potentially misrepresentative, proportion of potential EDB overwintering habitat that may have contributed to our inability to detect significant associations with EDB occupancy. This illustrates the shortcomings of using a monitoring study designed for a single species as a multispecies monitoring study and suggests that monitoring studies should be tailored toward individual species or those with similar ecologies and life histories or designed from the outset with multispecies inference in mind.
Eastern Indigo Snakes were more likely to occupy sandhills that had higher numbers of Gopher Tortoise burrows. This supports our initial hypothesis and corroborates the results of other studies noting the importance of Gopher Tortoise burrows for EIS winter retreat sites in the northern portion of their range (Stevenson et al., 2003, 2009; Hyslop et al., 2009). Additionally, our failure to detect EIS during our 3-yr study at 23 of our 40 sites, despite the presence of multiple tortoise burrows, strongly suggests that tortoise burrow presence alone is an insufficient predictor of suitable overwintering habitat. Our focus on only sandhills with tortoise burrows does not, however, directly assess the influence of tortoise burrow presence on EIS winter occupancy. Although this limits the inferential scope of our study, we have no reason to suspect that the sandhills we surveyed on protected lands do not comprise a representative sample of xeric sandhills supporting tortoise burrows throughout the Lower Altamaha River Watershed.
Eastern Indigo Snake occupancy also was significantly associated with the importance value of pine and amount of sandhill within 1 km. Pine importance had a strong negative correlation with oak importance (rs = −0.98), and sandhills with high oak importance often were fire-suppressed. Fire suppression may increase tree density/basal area and reduce ground cover abundance and diversity (Engstrom et al., 1984; Gilliam and Platt, 1999), thereby decreasing habitat quality for Gopher Tortoises. The low model support for the pine importance model, however, suggests the positive association between pine importance and EIS occupancy is conditional upon moderate–high numbers of tortoise burrows. Larger concentrations of sandhill habitat may correspond to greater amounts of potential overwintering habitat across the landscape that may support larger populations of EIS and increase connectivity among sandhills. The amount of sandhill was moderately correlated with number of tortoise burrows (rs = 0.54), however, and a post-hoc model with both of these covariates had equal ranking as either two single-variable models (max. ΔQAIC ≤ 0.39). Therefore, we were unable to fully assess the relative importance of these two covariates. Given the strong dependence of EIS on tortoise burrows (Stevenson et al., 2003, 2009; Hyslop et al., 2009), we suspect that tortoise burrow abundance has a more direct influence on EIS occupancy.
The strong model support of our site- and landscape-scale combination model suggests that factors operating at both scales influence winter EIS occupancy. Although receiving strong support and outperforming the null model, however, none of the landscape-scale covariates within this model were significant. The 95% CI for the amount of wetland and mixed forest within 1-km buffers had the least amount of overlap with zero, and their respective models were the highest ranked land-cover models, suggesting these two covariates were the most influential landscape-scale covariates within this model. The beta estimates for these two parameters suggests a negative association with wetlands and a positive association with mixed forest. The lack of association with land covers, however, particularly a negative association with wetlands, is not consistent with a previous study of EIS year-round habitat selection in southern Georgia. Hyslop et al. (2014) found that EIS used wetland, evergreen, mixed, and sparse (i.e., clearcut) land covers at both the study area- and home-range-scale more than deciduous and road/urban land covers. Wetland use was associated with foraging as 65% of observed foraging events were in wetlands; however, their study examined individual-level, year-round habitat selection, whereas our study examined population-level habitat associations during the winter when EIS foraging is likely reduced. Mismatches between individual- and population-level habitat selection studies may reflect scale-dependent (Bowyer and Kie, 2006; Mayor et al., 2009) or seasonally variable (Boyce et al., 2003) patterns of selection. Some researchers have reported different patterns of habitat selection across scales in individual-based studies of snake habitat selection (Harvey and Weatherhead, 2006; Lagory et al., 2009). Therefore, we suggest that differences in land cover associations between our study and Hyslop et al. (2014) are not contradictory per se but reflect habitat associations at different spatiotemporal scales, highlighting a greater need for understanding population-level EIS habitat associations.
There are multiple explanations for why we did not observe stronger relationships between EIS or EDB occupancy and our covariates. Our covariates may have had insufficient variability to allow us to detect an effect (Hartel et al., 2010; Fortin et al., 2012). All of our sites were located on protected lands which, although not necessarily providing suitable habitat, likely resulted in a greater degree of habitat quality both on and around our sites than might have been expected by randomly sampling sandhills from our study area. We do not consider this sufficient to completely explain our results, however, particularly for EIS. Habitat quality within our sites was highly variable, particularly with regard to tortoise burrow abundance and vegetation structure (Appendix 6). Additionally, occupancy is a relatively coarse-scale measurement of habitat suitability that does not incorporate variation in population size or probability of persistence. Habitat associations may appear stronger when compared to other population metrics such as abundance or fecundity. Finally, current patterns of EIS and EDB occurrence may reflect historical, rather than current, site and landscape conditions (Lindborg and Eriksson, 2004; Piha et al., 2007; Waldron et al., 2008; Halstead et al., 2014). Some species may show a lag in their response to anthropogenic habitat or landscape changes (Tilman et al., 1994; Metzger et al., 2009), particularly long-lived species that may persist for multiple generations following alterations to their habitat or the surrounding landscape (Kuussaari et al., 2009). Although our study was not designed to test this hypothesis, EIS and EDB may continue to use particular sandhills for overwintering, despite alterations to the sandhill or surrounding landscape. Alternatively, past land uses may have resulted in local extirpation but were followed by a degree of habitat recovery such that the current conditions appear suitable, but recolonization has not yet occurred.
Our study is the first to report estimates of EIS detection rates. Previous snake studies using visual encounter surveys have reported a wide range of detection rates (<0.10 to 0.70; Kery, 2002; Harvey, 2005; Christy et al., 2010; Sewell et al., 2012). Waldron et al. (2013b) used a winter survey methodology for EDB in coastal South Carolina similar to ours and reported a detection rate of 0.19 ± 0.11. In contrast, Steen et al. (2012b) reported lower monthly detection rates (≤0.17) for upland snakes in the southeastern Coastal Plain using box trapping from spring–fall, and their detection rates were ≤0.07 for EDB. We suspect that our EDB detection rates (0.22) were higher than those of Steen et al. (2012b), because we conducted our surveys when our study species were closely associated with visible habitat features (i.e., Gopher Tortoise burrows). Detection rates for EIS were higher than those for EDB, and we suggest potential explanations for this pattern. First, our survey seasons broadly overlapped the EIS breeding season that includes frequent interburrow mate-searching movements by adult males, male–male combat, and courtship (Speake et al., 1978; Stevenson et al., 2009; Hyslop et al., 2014). Gravid females also may increase their time spent basking in the late winter and early spring to raise their body temperature and facilitate vitellogenesis and/or egg development prior to oviposition (Blouin-Demers and Weatherhead, 2001; Shine, 2006; Lourdais et al., 2008). In contrast, EDB breeding activity in our study area is likely concluded by December (Timmerman and Martin, 2003; Hoss et al., 2011; Waldron et al., 2013a). Second, 35% of EDB detections were of snakes on the surface, whereas 85% of detections of EIS (snakes only) were on the surface. Third, EDB were rarely detected by shed skins, in contrast to EIS whose shed skins are regularly found near tortoise burrows (Diemer and Speake, 1983; Stevenson et al., 2003; Hyslop et al., 2009). Fourth, EDB are sit-and-wait ambush predators relying on cryptic coloration to avoid human interactions (Timmerman and Martin, 2003; Waldron et al., 2013b), whereas EIS are an actively foraging species and their semireflective black coloration increases their visibility. Finally, differences in detection may reflect differences in abundances between our study species, although we were unable to evaluate this hypothesis.
Our study species also differed in the factors affecting their detection. Air temperature had the strongest influence on EIS detection, consistent with their tendency for winter surface activity. Winter surface activity by EIS can occur under a wide range of weather conditions, including both clear and overcast skies with air temperatures ranging from 10°C to >25°C (Speake et al., 1978; Stevenson et al., 2003); yet our results indicate that detection is highest at warmer temperatures. Air temperature had comparatively little influence on EDB detection, perhaps because the majority of our detections were within tortoise burrows. Waldron et al. (2013b) also failed to detect an association between winter EDB detection rates and air temperature, although simulated data from radiotelemetered EDB showed a strong positive effect of air temperature on detection rate. In our study, survey date had the strongest influence on EDB detection rate, consistent with its late summer–fall breeding season and generally lower tendency for winter surface activity (Timmerman and Martin, 2003; Hoss et al., 2011; Waldron et al., 2013a). Both species showed evidence that our survey methodology influenced detection. Detecting an EIS at a site decreased the probability of detecting an EIS during a subsequent survey within a winter from 0.60 (95% CI = 0.36–0.80) to 0.33 (95% CI = 0.21–0.48). Detection rate for EDB varied from 0.10–0.28 across our three field observers. Although neither of these effects were significant, they highlight the importance of carefully designing survey methodologies to minimize heterogeneity in detection rates. Our results suggest that our survey methodology is better suited for detecting EIS and that additional survey effort or different survey dates may be required to substantially increase one's ability to detect EDB.
As anthropogenic demands on southeastern Coastal Plain landscapes continue, protected lands are likely to play an increasingly important role in species conservation. Care must be taken to ensure that protected lands provide the necessary ecological requirements to maintain species of conservation concern. Our results corroborate the previously noted importance of Gopher Tortoise burrows in xeric sandhills for EIS overwintering habitat, and highlight the importance of maintaining and restoring Gopher Tortoise populations. Management actions designed to benefit tortoise populations (e.g., prescribed burning) likely will benefit other longleaf pine ecosystem species as well (Means, 2006; Steen et al., 2013). Although our study was less well suited for examining the importance of these habitats for EDB overwintering habitat, our results do show that EDB will use xeric sandhills as overwintering habitat. Our study also illustrates that our winter survey methodology is effective in detecting EIS, and to a lesser extent EDB, and could be used to design future winter surveys for either species.
Acknowledgments
We thank the Georgia Department of Natural Resources for providing funding for this project and their personnel who helped make this study possible, particularly M. Harris, M. Elliott, J. Jensen, D. Mixon, and K. Morris. Riverbanks Zoo and E. Diebold provided additional funding through a Riverbanks Conservation Support Fund Grant. Additional funding was provided by The Orianne Society. A. Day and K. Briggs assisted in the field and H. Hall, M. Hodges, C. Martin, A. McGee, F. Snow, and W. Taylor assisted with various aspects of this study. L. Bailey provided advice on study design. Forest Investment Associates, Westervelt Company, and many private landowners provided us access to their properties. A. Allyn and K. McGarigal provided valuable suggestions for addressing spatial autocorrelation. This manuscript benefited greatly from the comments of T. Doan, P. Bartelt, and multiple anonymous reviewers. This research was conducted under United States Fish and Wildlife Service permit TE28025A-1 and Georgia Department of Natural Resources Scientific Collecting permit 29-WJH-13-143.
Literature Cited
Appendices
Appendix 1
Analysis of Residual Spatial Autocorrelation.—We tested for residual spatial autocorrelation (rSAC) by plotting correlograms of Moran's I calculated on our model residuals (Moore and Swihart, 2005). We calculated residuals following Moore and Swihart (2005) and created correlograms using the function correlog from the package “ncf” (v1.1-5, Bjornstad, 2013). We used distance bins of 800 m, which was twice the distance between our closest sites (Moore and Swihart, 2005), up to 20,000 m, which included the maximum distance-to-nearest-neighbor for all sites. We calculated autocovariate terms using the autocov_dist function in the “spdep” package (v0.5-88, Bivand, 2015) using both inverse Euclidean distance and inverse Euclidean distance squared weighting and used a neighborhood scheme following Bardos et al. (2015). We varied neighborhood size from 800–7,200 m based on assessments of correlograms. We then iteratively added different autocovariate terms to our models and selected the weighting method and neighborhood size that visually produced the greatest reduction in rSAC. We then reran our analyses including autocovariate terms in models where warranted.
We also assessed the degree to which our sites were ecologically independent. We considered sites to be ecologically independent if they were separated by a distance greater than the diameter of an EIS or EDB winter home range assuming a circular home range. We used published estimates of EIS winter (15 December through 15 March) home-range size from Hyslop et al. (2014), because these dates corresponded closely to our survey season. Specifically, we took the median of the upper 95% CI reported by sex and year (9.16 ha). This corresponded to a diameter of 342 m. Data on seasonal EDB home-range sizes are currently unavailable, but, given their low levels of winter surface activity (Waldron et al., 2013b), we assumed their winter home ranges would not exceed those of EIS. We plotted the frequency distribution of distances between nearest neighbors and calculated the proportion of nearest-neighbor distances less than or equal to the maximum winter home-range diameter. We calculated nearest-neighbor distance using both the distance between site centroids and the distance between site edges (i.e., sandhill polygon edges).
Appendix 2.
Moran's I correlograms of residuals from Eastern Indigo Snake (Drymarchon couperi, EIS) and Eastern Diamondbacked Rattlesnake (Crotalus adamanteus, EDB) single-season occupancy models. Autocovariates calculated using inverse Euclidean distance weighting were included in both EIS models (800 and 2,000 m neighborhood, respectively), whereas no autocovariates were included in both EDB models. Solid points represent significant (P < 0.05) residual spatial autocorrelation at a given 800 m distance bin. Correlograms depicted here are representative from the best-supported (ΔAIC/QAIC ≤ 2) models from each species: (A) EIS “global” site- and landscape-scale model containing number of tortoise burrows; (B) EIS model containing number of tortoise burrows; (C) EDB model containing the amount of evergreen forest in the 0.25 km buffer; and (D) EDB model containing importance value of pine.
Appendix 3.
Frequency distribution of distances-between-nearest-neighbors across our 40 survey sites in relation to the diameter of winter (15 December through 15 March) Eastern Indigo Snake (Drymarchon couperi) home ranges (median upper 95% CI = 342 m) from Hyslop et al. (2014). Nearest-neighbor distance was measured using site centroids (A) and the edge of each site's sandhill GIS polygon (B). Text boxes report the proportion of nearest-neighbor distances that were ≤ 342 m. Sites with nearest-neighbor distances > 342 m were considered ecologically independent.
Appendix 4.
Model selection results for the Eastern Indigo Snake (Drymarchon couperi) scale-selection analysis. Air temperature was used in the detection term for all models. An additive effect of year and a spatial autocovariate term were used in the occupancy terms for all models. Models are ranked according to their QAIC (AIC adjusted for overdispersion, c-hat = 1.61). Deviance (Dev) is calculated as −2*quasi-log-likelihood, K represents the number of parameters in the model, and wi is model weight.
Appendix 5.
Model selection results for the Eastern Diamondbacked Rattlesnake (Crotalus adamanteus) scale-selection analysis. A linear effect of survey date was used in the detection term for all models and an additive effect of year was included in the occupancy terms for all models. Models are ranked according to their AIC (c-hat = 1). Deviance (Dev) is calculated as −2*-log-likelihood, K represents the number of parameters in the model, and wi is model weight.
Appendix 6.
Summary statistics for sandhill- and landscape-scale covariates: minimum, 25th percentile, median, mean, 75th percentile, and maximum. All landscape-scale covariate names are reported as the land cover type followed by the buffer radius in km. All landscape-scale covariate values are reported as ha except patches (number of habitat patches) and edge density (density of habitat edges within each buffer).
