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.

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.

Results

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 (w_i = 0.82). The model with an effect of year (ΔAIC_C = 4.3, w_i = 0.10) and an additive model of the effects of phase and year (ΔAIC_C = 4.7, w_i = 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 (ΔAIC_C = 2.4). The inclusion of colonization and extinction probabilities in the latter model, however, meant it was penalized approximately four AIC_C 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 ΔAIC_C values of all models in the set were < 4.80 (Table 4), and four models had lower AIC_C 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. AIC_C, Akaike's Information Criterion adjusted for small sample size; Δ_i, the AIC_C of the model minus the AIC_C of the highest-ranked model; w_i, 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.

Discussion

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