Open Access
How to translate text using browser tools
21 April 2017 Relations between Environmental Attributes and Contemporary Occupancy of Threatened Giant Gartersnakes (Thamnophis gigas)
Eric C. Hansen, Rick D. Scherer, Erica Fleishman, Brett G. Dickson, Dave Krolick
Author Affiliations +

Habitat loss and modification are causing declines in the abundance and distribution of plant and animal species, yet robust information on which to base management and regulatory decisions for these species frequently is not available. Thamnophis gigas (Giant Gartersnakes), a species listed as threatened under the U.S. and California Endangered Species Acts, is strongly associated with aquatic ecosystems in the Great Central Valley (California), yet many aspects of its ecology remain poorly understood. We evaluated relationships between environmental attributes and occupancy of T. gigas and predicted the species' occupancy across ∼300,000 ha in the northern Central Valley. We trapped T. gigas at 159 sites and characterized land cover, land use, and soil type at each site. Occupancy of T. gigas was strongly and negatively associated with elevation and strongly and positively associated with canal density and the proportion of rice and perennial wetland. We also found a strong and previously undescribed association between occupancy and soil order. Estimated occupancy was over five times greater at sites underlain by alfisols, molisols, and vertisols than at sites underlain by entisols and inceptisols. We used the statistical associations between environmental variables and occupancy to predict occupancy at a spatial resolution and extent consistent with management of both T. gigas and regional land and water use.

Loss and modification of habitat are among the greatest threats to all taxonomic groups. Changes in land use and land cover in California's Great Central Valley (henceforth, Central Valley) since the mid-1800s exemplify such loss and modification in the world's major agricultural regions. To facilitate crop production, the Central Valley's once-expansive marshes and other wetlands were drained and channelized by reclamation levees, pumps, and canals (USFWS, 1999). An estimated 91% reduction in the area of California's wetlands has occurred since the 1780s, attributable to expansion of agriculture and urban areas (Dahl, 1990), with approximately 43% of freshwater wetlands in the Central Valley lost or converted since 1939 (Frayer et al., 1989). More recently, areas near Sacramento and other major cities are being converted from agriculture to urban and suburban housing. These changes have contributed to declines in the status of many wetland-associated animals and plants, yet robust ecological information on which to base management and regulatory decisions is not available for most species (Alford and Richards, 1999; Gibbon et al., 2000; Brooks et al., 2002). Agricultural and urban demands for water and other resources that may compete with species needs almost certainly will intensify.

Worldwide, reptiles are among the taxonomic groups for which ecological information is particularly sparse (Metrick and Weitzman, 1996; Gardner et al., 2007; Walsh et al., 2012). In the Sacramento Valley (northern Central Valley), lack of information on population dynamics and aspects of habitat quality for Thamnophis gigas (Giant Gartersnakes), a taxon listed as threatened under the U.S. and California Endangered Species Acts, impedes land-use decisions by diverse resource-management organizations. Thamnophis gigas is precinctive to the Central Valley and is more closely associated with aquatic environments than other gartersnakes in California (Fitch, 1940). The species inhabits low-gradient streams, valley floor wetlands, and marshes. It requires wetlands for foraging, upland areas for basking, upland burrows as summer shelter, and higher-elevation uplands for winter brumation (Hansen and Brode, 1980; USFWS, 1993, 1999; Hansen, 1998). All of these features were components of historical wetlands in the Central Valley; however, the distribution of this putative habitat for T. gigas is more extensive and more contiguous than the distribution of the species, suggesting that gradients of habitat quality or mechanisms of habitat selection may not be sufficiently understood.

Previous analyses identified environmental factors associated with the spatial distribution of T. gigas in the Central Valley and led to the generation of maps of habitat quality and predictions of occupancy for the species. On the basis of presence-only data, Halstead et al. (2010) reported that the likelihood of occurrence of T. gigas increased as proximity to rice fields and wetlands increased and as stream density decreased. In comparison, the results of models based on presence-only data tend to be less accurate than the results of models based on data on both presences and absences (Brotons et al., 2004). On the basis of presence and absence data from 24 sites, Halstead et al. (2014) concluded that patterns of T. gigas occupancy were best explained by distance to putative historic habitat (i.e., tule [Schoenoplectus acutus] marsh). Because T. gigas is absent from many areas of apparent habitat, Halstead et al. (2014) hypothesized that limitations to dispersal from historic habitat are constraining the species' contemporary spatial distribution. We agree that barriers to dispersal may be a factor, particularly over relatively large spatial extents; however, other factors may drive environmental heterogeneity that is associated with the distribution of the species at smaller extents.

The objective of this study was to evaluate hypothesized associations between the probability that a water body is occupied by T. gigas and attributes of the water body and adjacent lands. We used data from a much larger number of sites than did previous analyses, which gave us greater ability to evaluate potential causes of the patchy distribution of the species. We evaluated whether distance to historic tule marsh was associated with occupancy and assessed the strength of support for other hypotheses about components of habitat quality and selection for T. gigas. We used the statistical relations from our analyses to predict occupancy of T. gigas across a large portion of the northern Central Valley at a spatial extent consistent with regional management of the species and agricultural and urban expansion and operations.

Materials and Methods

Field Methods.—From 2009 through 2012, we conducted surveys of T. gigas in canals, sloughs, and wetlands in the Sacramento Valley (Fig. 1). We sampled 159 sites in the American and Yolo Basins and the southern Sutter Basin (Fig. 1). The study area covers at least 50% of the current range of T. gigas in the Sacramento Valley. Data for this analysis were derived from multiple studies, some of which included a random selection of sites and some of which reflected opportunistic sampling driven by perceived likelihood of snake occurrence. In all cases, we endeavored to maximize the range of representative landscape characteristics and distances among sites. The mean distance between sites was 15.54 km (range 0.07–44.14 km). At each site, consistent with methods described by Halstead et al. (2009, 2011), we sampled T. gigas along a line of 50 funnel traps spaced 10 m apart (total length of 500 m). We sampled public and, where landowners agreed, privately owned canals, ditches, wetlands, and drains in which water levels were sufficiently high to ensure that traps were continuously wet, thereby reducing the risk of desiccation or thermal stress for captured snakes. The traps were galvanized 4-mesh eel pots (Tackle Factory [Cuba Specialty Manufacturing], Fillmore, NY) modified to float following the procedures in Casazza et al. (2000). We placed traps along the edges of channels, streams, and associated vegetation, along which T. gigas tend to move and forage. In each of the 4 yr of the study, we sampled sites during two phases. The first phase was from early May through mid-June, and the second phase was from mid-July through mid-September. We chose these phases to encompass differences in life histories among sexes and age classes. During the first phase, males recently had left hibernacula and were exposed to traps while searching for mates and food (Coates et al., 2009). During the second phase, females had given birth and were exposed to traps while searching for food. We deployed traps for a minimum of 2 weeks at each site in a phase and maintained and checked all deployed traps daily.

Fig. 1.

Location of the Central Valley within California (top left), location of the study area relative to Sacramento (top right), and the sampling sites (bottom).


Environmental Data.—We derived nine environmental variables, eight continuous (Table 1) and one categorical, at each site. To derive four variables related to the land cover adjacent to sampled wetlands, we used ArcGIS (ArcGIS Desktop, release 10, Environmental Systems Research Institute, Redlands, CA) to delineate a 1,600-m buffer (the approximate distance that T. gigas move during their active season; E. Hansen, unpubl. data; USFWS, 1999) around each trap line and computed the proportion of 30-m cells within the buffer that fell within one of four classes: urban, rice or open water, row crops, and grassland (USDA, 2013). The spatial distributions of these classes changed little among the years of the study. Because all roads were classified by the U.S. Department of Agriculture (2013) as urban, the proportion of urban cover could be overestimated. Therefore, we retained clusters of urban cells and removed narrow strips of urban cells that usually represented roads (Theobald, 2013). We computed the density (length per unit area) of roads and canals within the 1,600-m buffer with the Line Density tool in ArcGIS. We then computed the mean density for all cells within the 1,600-m buffer. Elevation (which we measured at the centroid of the trap line) captured environmental attributes not represented by the other variables. We computed the distance to historic tule marsh with a digitized version of the map of Küchler (1977). We delineated a 10-m buffer around each transect, converted the buffer to a polygon, and used the Near tool in ArcGIS to compute the shortest distance between each transect and the polygon that represented the historic distribution of tule marsh. We hypothesized that occupancy is negatively associated with elevation, road density, distance to historic tule marsh, and the proportions of urban areas, grasslands, and row crops (Table 1). We also hypothesized that occupancy is positively associated with canal density and the proportions of rice and open water (Table 1; Halstead et al., 2010).

Table 1.

The continuous covariates we used in submodels of occupancy of Thamnophis gigas (Giant Gartersnakes) in the Sacramento Valley and sources of spatial data.


We used data from the National Cooperative Soil Survey (SSURGO) to classify the soils underlying each site as vertisols, mollisols, entisols-inceptisols, or alfisols. Entisols and inceptisols are separate orders in the SSURGO data, but we pooled the orders because the number of sites underlain by soils in either order was small and the two orders have similar topographic associations (e.g., both are young and are common in steep, rocky areas associated with rapid drainage). The majority of sites (N = 106) were located on vertisols; 28 and 11 sites were located on alfisols and mollisols, respectively, and 14 sites were on entisols and inceptisols. We hypothesized that historic hydrologic and geomorphologic processes drove soil formation in the valley, and that wetlands and their underlying soil types develop distinct hydrogeochemical signatures that affect habitat selection by T. gigas. These processes likely also affected the location of populations prior to conversion of wetlands to agriculture.

Statistical Analysis.—T. gigas, like many species of snakes, are difficult to capture or detect attributable to their wariness, periods of inactivity, and cryptic coloration (Lind et al., 2005; Durso et al. 2011; Halstead et al., 2013). Because the data were collected over 4 yr, a multiple-seasons occupancy model appeared to be appropriate; however, the number of parameters that must be estimated in a multiple-seasons model is greater than in a single-season model, and many of those parameters were tangential to our primary hypotheses. Therefore, to reduce the complexity of the model and avoid imprecise estimates of parameters, we analyzed the data with a single-season occupancy model (MacKenzie et al., 2002, 2006). Use of occupancy models to estimate probabilities of detection requires multiple surveys at each site (MacKenzie et al., 2002, 2006). We treated surveys during the two phases in each of the four years of the study as replicates, which resulted in eight possible surveys per site. Over 70% of the sites (N = 115) were sampled at least twice, and 38% of the sites (N = 60) were sampled four or more times. We considered the species to be detected at a given site in a given phase if at least one individual was captured in any trap.

The occupancy model included two parameters: pij, the probability that the species was detected at site i during survey j, and Ψi, the probability that the species occupied site i during any survey from 2009 through 2012. Prior to analyses, we developed sets of hypotheses describing the temporal variation in detection probability and the spatial variation in occupancy and converted these hypotheses into mathematical models. We fit models to the data in two stages (Lebreton et al., 1992). First, we combined models of detection probability with an intercept-only or null model of occupancy. Models of detection probability included fixed effects of year or phase. In models that included a fixed effect of phase, we allowed estimates of detection probability to vary among phases but not among years. In models that included fixed effects of year, we did not allow estimates of detection probability for the two phases within a year to vary, but allowed estimates to vary among years. We also fit models with additive and interactive effects of phase and year. Because we did not sample the same sites in all years and phases and weather conditions varied among sampling periods, we expected detection probability to vary among phases and years. We used Akaike's information criterion adjusted for small sample size (AICC) and Akaike weights (wi) to evaluate the models of detection probability and in all other model selection procedures (Burnham and Anderson, 2002). We considered models with AICC values within 2 points of one another to have similar statistical support (Burnham and Anderson, 2002).

In the second stage, we combined all models of occupancy with the highest-ranked models of detection probability (the model with the lowest AICC value and any models within 2 AICC points). The models of occupancy included the nine environmental variables described above. Initially, we fit models to the data with an intercept-only model of occupancy or an intercept and a single environmental variable. We identified the variables for which the 95% confidence intervals (CIs) around the estimates of the regression coefficients did not include 0 and included those variables in models of occupancy with three variables. We limited the number of environmental variables to three to prevent inclusion of correlated variables in the same model, and we did not include two variables in the same model if their correlation coefficient was > |0.60|. Similarly, we did not include multiple land-cover variables in the same model because they are inherently correlated (i.e., an increase in the proportion of one land-cover type decreased the proportions of the other land-cover types). We evaluated the relative strengths of association of the environmental variables with occupancy on the basis of their occurrence in the highest-ranked models and the magnitudes of the estimates of their regression coefficients. We report the estimates of regression coefficients from the highest-ranked model in which each covariate was included and 95% CIs for all point estimates. We tested the goodness-of-fit of the highest-ranked model to the data with the bootstrap method of MacKenzie and Bailey (2004).

An assumption of the occupancy model is that the occupancy status of each sample unit does not change over the duration of sampling (the closure assumption). Because the sample units in this study were small relative to the movement distances of T. gigas (USFWS, 1999; E. Hansen, unpubl. data), snakes could have emigrated from or immigrated to sites between surveys. In addition, the occupancy status of a given site may have changed during the 4 yr of sampling. Under some conditions, violations of the closure assumption can bias estimates of occupancy (MacKenzie et al., 2006). Therefore, in addition to assessing the goodness-of-fit of the model, we organized the data as though they were collected under Pollock's robust design (Pollock, 1982) and used the method of Rota et al. (2009) to test whether the assumption of closure was met. We treated the years as seasons and fit the multiple-seasons occupancy model (MacKenzie et al., 2003, 2006) to the data. The multiple-seasons occupancy model allows for changes in site-level occupancy status among seasons and includes two parameters that represent the processes that lead to changes in status: γy, the probability that a site is colonized between seasons y and y + 1; and ϵy, the probability that the species becomes extinct at a site between seasons y and y + 1. We fit three versions of the multiple-seasons model to the reorganized data. In the first model, estimates of γy and ϵy varied among years. In the second model, γy and ϵy were constant among years. In the third model, γy and ϵy were fixed to 0 (to represent the closure assumption).

We standardized all of the continuous environmental variables prior to fitting occupancy models to the data; however, because one of our objectives was to predict occupancy of T. gigas at locations that we did not sample, we did not standardize the variables with the set of values from the 159 sites. Instead, we standardized the variables and predicted occupancy across a larger area (Dickson et al., 2013) delineated by the outer boundary of the set of 1,600-m buffers around each site. We placed a 1,600-m buffer around each 30-m cell within the boundary and derived values of the eight environmental variables for each cell.

Predictions of Occupancy.—We used model-averaged estimates of the slope and regression coefficients from the models with three environmental variables to predict occupancy of T. gigas at unsampled locations. We assigned a value of 0 to regression coefficients for environmental variables not included in models (Burnham and Anderson, 2002). The range of values of environmental variables in the area in which we predicted habitat quality is broader than the range of values among the sites that we surveyed (Table 2). Therefore, some predictions were based on values of variables that were outside the range of values in the field data.

Table 2.

The range of values of continuous covariates across all sites and across the area in which we predicted occupancy for Thamnophis gigas (Giant Gartersnakes) in the Sacramento Valley.


To conduct a preliminary test of the spatially explicit predictions, we used data from the California Natural Diversity Database (CNDDB, 2015). The data, which are presence-only, include coordinates of observations of T. gigas from across the Central Valley. Although most of the records from the CNDDB are from the last 20–30 yr, some records are older, and T. gigas may have been extirpated from areas in which they were observed more than three decades ago. Despite these shortcomings, however, the data are a useful and independent means of evaluating the predictions of our models. We used ArcGIS to extract the value of habitat quality from each cell in which T. gigas was recorded and calculated the percentage of values that exceeded 0.50.


We detected T. gigas during at least one survey at 82 sites (52%). The highest-ranked model of detection probability included an interaction between the effects of phase and year and had nearly all the model weight (wi = 0.82). The model with an effect of year (ΔAICC = 4.3, wi = 0.10) and an additive model of the effects of phase and year (ΔAICC = 4.7, wi = 0.08) on detection probability had considerably less support in the data, and we did not include them in the second stage of model fitting. Model-averaged estimates of detection probability for a given phase and year ranged from 0.10 (CI 0.01–0.47) to 0.81 (0.61–0.92). Six of the eight estimates were ≥ 0.30, and three of the estimates were > 0.70 (Table 3). The highest-ranked model fit the data well (ĉ = 0.85).

Table 3.

Model-averaged estimates of detection probability for both phases in each of 4 yr, and their lower and upper 95% confidence intervals, for Thamnophis gigas (Giant Gartersnakes) in the Sacramento Valley. Phase 1 was from early May through mid-June, and phase 2 was from mid-July through mid-September.


Our results indicated that the assumption of closure was met. The highest-ranked multiple-seasons model included colonization and extinction probabilities that were fixed to 0. The model in which colonization and extinction probabilities were estimated and constrained to be equal among years was ranked second and also was supported by the data (ΔAICC = 2.4). The inclusion of colonization and extinction probabilities in the latter model, however, meant it was penalized approximately four AICC units. Because the deviance of the model accounted for less than half the penalty, inclusion of extinction and colonization probabilities may not be important (Anderson, 2008). Furthermore, estimates of extinction probability from the latter model were low (ε̂y = 0.12 [0.03–0.37]), and estimates of colonization probability were < 0.01, although highly imprecise. Of the models of occupancy that included a single environmental variable, the highest ranked included canal density. The association between occupancy and canal density was positive. The estimated 95% CI for the regression coefficient of canal density did not include 0. Estimates of 95% CIs for four additional environmental variables also did not include 0, and we included these variables in three-variable models of occupancy. The latter variables were elevation (negative association with occupancy), soil order (lower occupancy for entisols and inceptisols), proportion of row crops (negative), and proportion of rice or open water (positive). There was no support for the model that included distance to historic tule marsh, and the 95% CI around the estimate of the regression coefficient from the model overlapped 0.

The strengths of support for the three-variable models of occupancy were similar. The ΔAICC values of all models in the set were < 4.80 (Table 4), and four models had lower AICC values than the highest-ranked single-variable model. Elevation and canal density were the most strongly supported variables (Table 4). Of the four highest-ranked models, they were in four and three of the models (Table 4). Associations of occupancy with proportion of rice or open water, proportion of row crops, and soil order had weaker support, although 95% CIs around estimates of their regression coefficients did not include 0 or slightly overlapped 0 (Tables 4, 5). Estimates of regression coefficients from the models that included soil order indicated that occupancy was lower at sites on entisols or inceptisols than at sites on the other three soil orders. Estimated occupancy in areas underlain by entisols and inceptisols was more than five times lower than estimated occupancy in areas on alfisols, mollisols, and vertisols.

Table 4.

Model selection results for the three-covariate models of occupancy of Thamnophis gigas (Giant Gartersnakes) in the Sacramento Valley. AICC, Akaike's Information Criterion adjusted for small sample size; Δi, the AICC of the model minus the AICC of the highest-ranked model; wi, Akaike weight; k, number of parameters; −2*log(L), −2 times the value of the likelihood at its maximum; Ψ(.) , intercept-only model. Environmental covariates are defined in Table 1.


Table 5.

Relative strengths of association between occupancy of Thamnophis gigas (Giant Gartersnakes) and each covariate included in the four highest-ranked, three-covariate models of occupancy. Estimates of regression coefficients from the highest-ranked model with each covariate and model-averaged estimates of regression coefficients are also provided. Model-averaged estimates were used in the map of predicted occupancy.


Predicted probabilities of occupancy varied from 0.00 to 0.94, and the highest values were concentrated in the Colusa, Sutter, American, and Yolo Basins (Fig. 2). These basins were separated by areas with low values of predicted occupancy that largely reflect the entisols and inceptisols adjacent to the Sacramento River (Figs. 2, 3). Over 70% of the presences of T. gigas from the CNDDB occurred in cells with predicted occupancy > 0.50 (Figs. 2, 3).

Fig. 2.

Predicted probability of Thamnophis gigas (Giant Gartersnakes) occupancy in the northern Central Valley of California. Blue, lowest predicted occupancy; red, highest predicted occupancy; black dots, presences of T. gigas from the California Natural Diversity Database. AB, American Basin; CB, Colusa Basin; NB, Natomas Basin (southern American Basin); SB, Sutter Basin; YB, Yolo Basin.


Fig. 3.

Occupancy of Thamnophis gigas (Giant Gartersnakes) in the Natomas Basin. Transparent purple and orange overlays represent the Natomas and American Basins, respectively, and the blue dots represent presences of T. gigas from the California Natural Diversity Database.



Thamnophis gigas use canals as refugia in spring before rice fields have been inundated and in autumn as rice fields are being dewatered after the harvest (Hansen et al., 2015). Canals also may be necessary components of T. gigas habitat given the dynamic agricultural practices in the Sacramento Valley, where rice fields may be fallowed or converted to other crops as market conditions change. Rice fields also may be fallowed if the value of water increases and farmers choose to sell their water rights. Although T. gigas populations are not likely to persist over long periods without access to inundated rice fields, canals may provide sufficient resources during relatively short periods of fallowing. Rice fields are believed to function as surrogates for the wetlands that formerly were core habitat for the species (Hansen, 1998; Halstead et al., 2014) and to serve as essential components of habitat during the active season of T. gigas. The rice fields' shallow, warm waters produce high concentrations of aquatic prey, and the rice plants provide cover from predators (Hansen, 1998; Halstead et al., 2010). Nevertheless, occupancy of rice fields is heterogeneous, even in areas where no historical or contemporary barriers to movement are apparent.

Many mechanisms may explain the negative association between occupancy of T. gigas and elevation, including the development of different prey assemblages associated with rapid upland drainage and lowland perennial marsh, the relative lack of wetlands and other aquatic systems on the hillsides adjacent to the valley floor, and the increased energetic cost of dispersal from low elevations to wetlands at higher elevations. In addition, different soil orders are more common on the valley floor and on the surrounding hillsides. Entisols and inceptisols are recently formed, coarse-loamy soils with a geomorphic history of rapid drainage associated with elevation change and coarse sediment deposition by floodwaters. By contrast, alfisols, mollisols, and vertisols generally occur in relatively flat areas (e.g., steppes and prairies). These three clay-rich, fertile, and fine-textured soil orders are associated with hydric conditions that result from frequent flooding and persistent inundation. They also occur beneath historic tule marsh, and support the majority of rice agriculture on the floor of the Sacramento Valley.

Therefore, effects of elevation may be confounded with effects of soil order. Although the effect of elevation was more strongly supported than soil order in the highest ranked model in the original model set, the greater number of parameters in the original model of soil order may have been the cause. In the original model, we treated soil as a fixed effect with four orders, which meant the model had two more parameters than the models with elevation or other continuous covariates. Model-selection metrics such as AICC generally penalize models with larger numbers of parameters (Burnham and Anderson, 2002). A posteriori, we fit a second version of the soil model in which occupancy at sites underlain by vertisols, alfisols, and mollisols was equal but differed from occupancy at sites underlain by entisols and inceptisols. We fit this second model because estimates of regression coefficients from the initial model indicated that occupancy did not differ among sites on vertisols, alfisols, and mollisols. The revised model also had the same number of parameters as the models with continuous covariates. The revised model became the highest ranked, and the estimate of the regression coefficient indicated much lower occupancy at sites underlain by entisols and inceptisols. These results are consistent with the hypothesis that underlying soil order is associated with occupancy of T. gigas, and they provide one explanation for observed, finer-resolution patterns of the species' distribution in the study area. For instance, our results suggest that occupancy of T. gigas differs in areas underlain by coarse, permeable alluvium along the Sacramento River through the southwestern and southern portions of the Natomas Basin (near the Sacramento International Airport) and in the eastern portion of the basin (Fig. 3). Historic and contemporary observations of T. gigas in the southwestern and southern portions of the basic are rare, whereas there are many observations of the species elsewhere in the basin.

The association between soil order and occupancy by T. gigas may be explained by chemotactic selection. For example, different soil orders may be characterized by distinct geochemical attributes that T. gigas recognizes. Use of chemoreception is well documented in gartersnakes. For example, male gartersnakes use chemoreception to locate females that excrete pheromones through lipids in the skin (Ford and Low, 1984; Ford and O'Bleness, 1986; Mason, 1993; LeMaster et al., 2001) and to discriminate between mated and unmated females (O'Donnell et al., 2004). Pheromones also may be used by gartersnakes to locate traditional dens during autumnal migrations (Costanzo, 1989; Lemaster et al., 2001). Although areas underlain by entisols and incepticols often are leveled to maintain persistent surface water for use in rice agriculture, they still may retain the historic chemical signature of rapidly draining soils, and our results suggest that T. gigas may avoid these areas.

Previous analyses cited the estimated extent of historic tule marsh as the variable most strongly associated with T. gigas occupancy (Halstead et al., 2014), but the results of our analyses did not support that hypothesis. We used data from a greater number of sites than Halstead et al. (2014) (N = 159 vs. N = 24), and our sites were concentrated in a smaller area. In addition, ∼74% of our sites were within the area of predicted historic tule marsh. At the coarser resolution examined by Halstead et al. (2014), mechanisms related to dispersal appeared to be more relevant. They evaluated covariates that were measured at multiple resolutions, and models that included the distance to historic tule marsh were the most strongly supported. The higher density of sites in our study allowed us to evaluate mechanisms that operate at finer resolutions. At this resolution, canal density, the proportion of adjacent rice agriculture and wetlands, and underlying soils appear to be stronger drivers of T. gigas occupancy.

Legally mandated conservation actions for T. gigas are included within numerous regional land-use plans (e.g., City of Sacramento et al., 2003; Yolo County, 2015). Some plans include habitat area or configuration requirements (City of Sacramento et al., 2003; Yolo County, 2015) or simple assumptions about attributes of vegetation structure and composition that equate to habitat for the species (Yolo County, 2015). In many cases, these assumptions have not been validated by robust analyses of comprehensive empirical data and are not consistent with our results, suggesting a higher degree of heterogeneity in T. gigas habitat in the Central Valley.

The spatial scale of our results is highly relevant to regional resource management (e.g., City of Sacramento et al., 2003; USFWS 2015). For example, our results facilitate the identification of individual habitat patches within basins and barriers to migration of T. gigas among basins. Results from analyses of nuclear DNA suggested similar connections or barriers among basins (Wood et al., 2015). For example, both occupancy and genetic results indicated that T. gigas populations in the American Basin (including the Natomas Basin, or southern American Basin) are relatively well-connected to basins on the east side of the Sacramento River (e.g., Sutter Basin and Butte Basin) but less connected to populations on the west side of the Sacramento River (e.g., Yolo Basin) (Wood et al., 2015). Together or individually, such data on occupancy and population genetics can be used to estimate current and future patterns of movement, gene flow, and genetic differentiation given different scenarios of environmental change (e.g., McRae et al., 2008). This information also is useful for identifying the most effective spatial scales for taking conservation actions, such as mitigation or translocation.

Although we used data from multiple years, we fit a single-season occupancy model because our primary objective was to identify correlates of occupancy. We also were aware that the low detection probability of T. gigas may require pooling data across years to generate estimates that are reasonably precise. We tested the assumption of closure, and the results of our test supported use of the single-season model; however, we acknowledge that the low capture probabilities preclude strong inference regarding closure. Based on our experience with the species and the results of other studies, we do not think that violations of the closure assumption biased the inferences from this study. For example, the positive associations between occupancy and canal density and proportion of rice in our study were consistent with the results of Halstead et al. (2010), who collected data in the same general locations but used presence-only data and different analytical methods (factor analyses; Calenge and Basille, 2008).

The precision of our estimates was strengthened by the detection probabilities in 2009 and 2010 (Table 3), which are relatively high for snakes (Durso et al., 2011). During those years, we sampled wetlands in the core of the current distribution of T. gigas, where we expected densities to be moderate or high (e.g., USFWS, 1999, 2015; CNDDB, 2015). In 2011 and 2012, we sampled a higher proportion of wetlands outside the core, where we expected that densities of T. gigas would be much lower. Our expectations about the gradient of density appeared to be correct, and, consequently, estimates of detection probability declined and were generally less precise (Table 3). These results suggest that occupancy studies based on data from aquatic trapping may lead to imprecise estimates when T. gigas densities are low, and other sampling methods may be more effective. We are currently investigating the usefulness of other methods, and preliminary results suggest that environmental DNA (E. Hansen and G. Schumer, unpubl. data) and scent detection dogs (E. Hansen and H. T. Harvey and Associates, unpubl. data) may be feasible alternatives where the probability of detecting T. gigas is low.

Although our work increases the resolution of estimates of T. gigas occupancy, our predictions of occupancy were based in part on values of covariates outside the range of our field data. Additionally, there are gaps and regional inconsistencies in the SSURGO soil data, and the number of trap lines in areas underlain by entisols and inceptisols was small. Therefore, we suggest that our predictions be evaluated with additional data. We also suggest that future work emphasize identification of soil-chemistry metrics that can be measured in the field and that are indicative of biologically meaningful soil properties. Such information could facilitate rapid assessment of field conditions, improve the capacity to predict occupancy, and facilitate reliable, fine-resolution predictions of habitat quality throughout the species' range. Such knowledge might improve the ability to manage T. gigas at a regional level and across its range.

Acknowledgments.—We thank the dedicated field personnel that made data collection possible. Two anonymous reviewers provided comments that improved the quality of this manuscript. Dylan Harrison-Atlas provided additional assistance with spatial analyses. Support for this work was provided by Brookfield Natomas LLC, the Sacramento Area Flood Control Agency, the Natomas Basin Conservancy, and the Conaway Preservation Group. Reclamation District 1000, Reclamation District 1001, and Reclamation District 1500 provided invaluable assistance with site access. This study was conducted under U.S. Fish and Wildlife Service Recovery Permit 10(a) (1) (A) ESA TE-018177-7 and Department of Fish and Wildlife Scientific Collection Permit SC-003881.

Literature Cited


Alford, R. A., and S. J. Richards. 1999. Global amphibian declines: a problem in applied ecology. Annual Review of Ecology and Systematics 30:133–165. Google Scholar


Anderson, D. R. 2008. Model Based Inference in the Life Sciences: A Primer on Evidence. New York: Springer. Google Scholar


Brooks, T. M., R. A. Mittermeier, C. G. Mittermeier, G. A. B. da Fonseca, A. B. Rylands, W. R. Konstant, P. Flick, J. Pilgrim, S. Oldfield, G. Magin, and C. Hilton-Taylor. 2002. Habitat loss and extinction in the hotspots of biodiversity. Conservation Biology 16:909–923. Google Scholar


Brotons, L., W. Thuiller, M. B. Araújo, and A. H. Hirzel. 2004. Presence-absence versus presence-only modelling methods for predicting bird habitat suitability. Ecography 27:437–448. Google Scholar


Burnham, K. P., and D. R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. Springer-Verlag, USA. Google Scholar


Calenge, C., and M. Basille. 2008. A general framework for the statistical exploration of the ecological niche. Journal of Theoretical Biology 252:674–685. Google Scholar


CNDDB (California Natural Diversity Database). 2015. Thamnophis gigas. California Department of Fish and Wildlife. California, USA [Internet]. Available from: Accessed 9 March 2015. Google Scholar


Casazza, M. L., G. D. Wylie, and C. J. Gregory. 2000. A funnel trap modification for surface collection of aquatic amphibians and reptiles. Herpetological Review 31:91–92. Google Scholar


City of Sacramento, Sutter County, and Natomas Basin Conservancy. 2003. Final Natomas Basin Habitat Conservation Plan. California, USA [Internet]. Available from: Archived by WebCite at on 6 October 2016. Google Scholar


Coates, P. S., G. D. Wylie, B. J. Halstead, and M. L. Casazza. 2009. Using time-dependent models to investigate body condition and growth rate of the giant gartersnake. Journal of Zoology 279:285–293. Google Scholar


Costanzo, J. P. 1989. Conspecific scent trailing by garter snakes (Thamnophis sirtalis) during autumn: further evidence for the use of pheromones in den location. Journal of Chemical Ecology 15:2531–2538. Google Scholar


Dahl, T. E. 1990. Wetlands Losses in the United States 1780s to 1980s. U.S. Department of the Interior, U.S. Fish and Wildlife Service, USA. Google Scholar


Dickson, B. G., S. E. Sesnie, E. Fleishman, and D. S. Dobkin. 2013. Identification of habitat and assessment of habitat quality for conservation of terrestrial animals. Pp. 149–173in L. Craigheadand C. Convis (eds.), Conservation Planning: Shaping the Future. Esri Press, USA. Google Scholar


Durso, A. M., J. D. Willson, and C. T. Winne. 2011. Needles in haystacks: estimating detection probability and occupancy of rare and cryptic snakes. Biological Conservation 144:1508–1515. Google Scholar


Fitch, H. S. 1940. A biogeographical study of the ordinoides artenkreis of garter snakes (genus Thamnophis). University of California Publications in Zoology 44:1–150. Google Scholar


Ford, N. B., and J. R. Low. 1984. Sex pheromone source location by garter snakes: a mechanism for detection of direction in nonvolatile trails. Journal of Chemical Ecology 10:1193–1199. Google Scholar


Ford, N. B., and M. L. O'Bleness. 1986. Species and sexual specificity of pheromone trails of the garter snake, Thamnophis marcianus. Journal of Herpetology 20:259–262. Google Scholar


Frayer, W. E., D. D. Peters, and H. R. Pywell. 1989. Wetlands of the California Central Valley: Status and Trends. U.S. Fish and Wildlife Service, USA. Google Scholar


Gardner, T. A., J. Barlow, and C. A. Peres. 2007. Paradox, presumption and pitfalls in conservation biology: the importance of habitat change for amphibians and reptiles. Biological Conservation 138:166–179. Google Scholar


Gibbon, J. W., D. E. Scott, T. J. Ryan, K. A. Buhlmann, T. D. Tuberville, B. S. Metts, J. L. Greene, T. Mills, Y. Leiden, S. Poppy et al. 2000. The global decline of reptiles, déjà vu amphibians. BioScience 50:653–666. Google Scholar


Halstead, B. J., G.D. Wylie, P. S. Coates, and M. L. Casazza. 2009. The U.S. Geological Survey quantitative adaptive survey protocol for the giant gartersnake (Thamnophis gigas). U.S. Geological Survey, Western Ecological Research Center, USA. Google Scholar


Halstead, B. J., G. D. Wylie, and M. L. Casazza. 2010. Habitat suitability and conservation of the giant gartersnake (Thamnophis gigas) in the Sacramento Valley of California. Copeia 2010:591–599. Google Scholar


Halstead, B. J., G. D. Wylie, P. S. Coates, and M. L. Casazza. 2011. Bayesian adaptive survey protocols for resource management. Journal of Wildlife Management 75:450–457. Google Scholar


Halstead, B. J., G. D. Wylie, and M. L. Casazza. 2013. Efficacy of trap modifications for increasing capture rates of aquatic snakes in floating aquatic funnel traps. Herpetological Conservation and Biology 8:65–74. Google Scholar


Halstead, B. J., 2014. Ghost of habitat past: historic habitat affects the contemporary distribution of giant garter snakes in a modified landscape. Animal Conservation 17:144–153. Google Scholar


Hansen, G. E. 1998. Cherokee Canal sediment removal project post-construction giant garter snake (Thamnophis gigas) surveys. Final Report for California Department of Water Resources. Contract No. B-81535. Google Scholar


Hansen, G. E., and J. M. Brode. 1980. Status of the giant garter snake, Thamnophis couchii gigas. Inland Fisheries Endangered Species Program Special Publication 80-5:1–14. Google Scholar


Hansen, E. C., R. D. Scherer, G. C. White, B. G. Dickson, and E. Fleishman. 2015. Estimates of survival probability from two populations of giant gartersnakes in California's Great Central Valley. Copeia 103:1026–1036. Google Scholar


Küchler, A.W. 1977. Appendix: the map of the natural vegetation of California. Pp. 909–938in M. G. Barbourand J. Major (eds.), Terrestrial Vegetation of California. John Wiley and Sons, USA. Google Scholar


Lebreton, J. D., K. P. Burnham, J. Clobert, and D. R. Anderson. 1992. Modeling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecological Monographs 62:67–118. Google Scholar


Lemaster, M. P., I. T. Moore, and R. T. Mason. 2001. Conspecific trailing behaviour of red-sided garter snakes, Thamnophis sirtalis parietalis, in the natural environment. Animal Behaviour 61:827–833. Google Scholar


Lind, A. J., H. H. Welsh Jr., and D. A. Tallmon. 2005. Garter snake population dynamics from a 16-year study: considerations for ecological monitoring. Ecological Applications 15:294–303. Google Scholar


MacKenzie, D. I., and L. L. Bailey. 2004. Assessing the fit of site-occupancy models. Journal of Agricultural, Biological, and Environmental Statistics 9:300–318. Google Scholar


MacKenzie, D. I., J. D. Nichols, G. B. Lachman, S. Droege, J. A. Royle, and C. A. Langtimm. 2002. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83:2248–2255. Google Scholar


MacKenzie, D. I., J. D. Nichols, J. E. Hines, M. G. Knutson, and A. B. Franklin. 2003. Estimating site occupancy, colonization, and location extinction when a species is detected imperfectly. Ecology 84:2200–2207. Google Scholar


Mackenzie, D. I., J. D. Nichols, J. A. Royle, K. H. Pollock, L. L. Bailey, and J. E. Hines. 2006. Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence. Academic Press, USA. Google Scholar


Mason, R. T. 1993. Chemical ecology of the red-sided garter snake, Thamnophis sirtalis parietalis. Brain, Behavior and Evolution 41:261–268. Google Scholar


McRae, B. H., B. G. Dickson, T. H. Keitt, and V. B. Shah. 2008. Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology 89:2712–2724. Google Scholar


Metrick, A., and M. L. Weitzman. 1996. Patterns of behavior in endangered species preservation. Land Economics 72:1–16. Google Scholar


O'Donnell, R. P., N. B. Ford, R. Shine, and R. T. Mason. 2004. Male red-sided garter snakes, Thamnophis sirtalis parietalis, determine female mating status from pheromone trails. Animal Behaviour 68:677–683. Google Scholar


Pollock, K. H. 1982. A capture-recapture design robust to unequal probability of capture. Journal of Wildlife Management 46:752–757. Google Scholar


Rota, C. T., R. J. Fletcher Jr., R. M. Dorazio, and M. G. Betts. 2009. Occupancy estimation and the closure assumption. Journal of Applied Ecology 46:1173–1181. Google Scholar


Theobald, D. M. 2013. A general model to quantify ecological integrity for landscape assessments and U.S. application. Landscape Ecology 28:1859–1874. Google Scholar


USDA (U.S. Department of Agriculture). 2013. National Agricultural Statistics Service Cropland Data Layer [Internet]. Available from: Archived by WebCite at on 6 October 2016. Accessed 2013. Google Scholar


USFWS (U.S. Fish and Wildlife Service). 1993. Endangered and threatened wildlife and plants: determination of threatened status for the giant garter snake. Federal Register 58:54053–54066. Google Scholar


USFWS (U.S. Fish and Wildlife Service). 1999. Draft Recovery Plan for the Giant Garter Snake (Thamnophis gigas). U.S. Fish and Wildlife Service, USA. Google Scholar


USFWS (U.S. Fish and Wildlife Service). 2015Revised Draft Recovery Plan for the Giant Garter Snake (Thamnophis gigas). U.S. Fish and Wildlife Service, Pacific Southwest Region, USA. Google Scholar


Walsh, J. C., J. E. M. Watson, M. C. Bottrill, L. N. Joseph, and H. P. Possingham. 2012. Trends and biases in the listing and recovery planning for threatened species: an Australian case study. Oryx 47:134–143. Google Scholar


Wood, D. A., B. J. Halstead, M. L. Casazza, E. C. Hansen, G. D. Wylie, and A. G. Vandergast. 2015. Defining population structure and genetic signatures of decline in the giant gartersnake (Thamnophis gigas): implications for conserving threatened species within highly altered landscapes. Conservation Genetics 16:1025–1039. Google Scholar


Yolo County HCP / NCCP Joint Powers Authority. 2015. Second Administrative Draft Yolo Habitat Conservation Plan and Natural Community Conservation Plan [Internet]. Available at Archived by WebCite at on 5 October 2016. Google Scholar
Copyright 2017 Society for the Study of Amphibians and Reptiles
Eric C. Hansen, Rick D. Scherer, Erica Fleishman, Brett G. Dickson, and Dave Krolick "Relations between Environmental Attributes and Contemporary Occupancy of Threatened Giant Gartersnakes (Thamnophis gigas)," Journal of Herpetology 51(2), 274-283, (21 April 2017).
Accepted: 1 February 2017; Published: 21 April 2017
Back to Top