Skip to main content
Dryad logo

Variation in behavior drives multiscale responses to habitat conditions in timber rattlesnakes (Crotalus horridus)


Hoffman, Andrew et al. (2022), Variation in behavior drives multiscale responses to habitat conditions in timber rattlesnakes (Crotalus horridus), Dryad, Dataset,


Variations in both the behavior of wildlife and the scale at which the environment most influences the space use of wild animals (i.e., scale of effect) are critical, but often overlooked in habitat selection modeling. Ecologists have proposed that biological responses happening over longer time frames are influenced by environmental variables at larger spatial scales, but this has rarely been empirically tested. Here, we hypothesized that long-term patterns of behavior (i.e., lasting multiple weeks to months) would be associated with larger scales of effect than more sporadic behaviors. We predicted site use by 43 radio-telemetered timber rattlesnakes (Crotalus horridus) exhibiting four distinct, time-varying behaviors (foraging, digestion, ecdysis, and gestation) using remotely-sensed environmental variables related to forest structure and landscape topography. Among sites used by snakes, warmer temperatures and higher levels of forest disturbance were predictive of behaviors dependent on thermoregulation including gestation and ecdysis while more moderate temperatures and drier, more oak-dominated sites were predictive of foraging. Long-term behaviors were associated with larger spatial scales across most variables, supporting our hypothesis that the scale at which habitat selection occurs is linked to the temporal scale of relevant behaviors. Management recommendations based on single-scale models of habitat use that do not account for fine-scale variations in behavior may obscure the importance of potentially limiting habitat features needed for infrequent behaviors that are important for growth and reproduction of this and related species.


Telemetry – We captured individual timber rattlesnakes and surgically implanted intraperitoneal transmitters (Holohil SI-2T) and then tracked snakes between dawn and dusk, regularly shifting the order in which individuals were tracked daily to avoid systematic bias. Upon visually locating a snake, we recorded the location with a Global Positioning System (Garmin GPSmap 64s) to an estimated < 5m spatial accuracy. Raw snake location data in the form of X, Y coordinates are not included in this dataset as they are sensitive and not necessary to run the models used in our analysis. Instead, each unique point location is given an alpha numeric code in the "Point" column. Each row represents a unique snake observation for a given date (indicated by the "Date" column) and individual snake (indicated by the "Snake" column). The "dat_use_5m.csv" also contains paired random points that do not represent snake observations (indicated by a 0 in the "Used" column).

Behavior Classification – We identified four behavioral and physiological states (hereafter referred to as “behaviors”): foraging, digestion, ecdysis, and gestation. We considered snakes to be foraging when they were observed in a stereotyped foraging posture. Our observations of snakes digesting a recent meal are limited by the difficulty in detecting a small bolus in a large-bodied snake, but we noted this whenever possible. These observations are, therefore, inherently biased toward snakes digesting larger, more obvious meals and may be prone to a high false negative detection rate. We noted observations of snakes in ecdysis as indicated by snakes having blue-gray eyes and dusky skin coloration. We confirmed observations of gestating females by sonogram shortly after snake emergence in April and May. We excluded all locations in which a snake was not visible and presumably underground or actively moving. When not moving and not clearly participating in the aforementioned behaviors, we designated snake behavior as resting. We excluded snake observations from analyses when there was uncertainty about behavioral classification. For each snake observation, we indicate which of these four behaviors, if any, were observed in their respective behavioral columns ("Forage", "Ecdysis", "Digest", and "Gestate"; a 1 indicates the behavior was observed). The dat_use_5m.csv file does not contain any behavioral columns as it is only to be used when running the general site use model.

Geospatial Habitat Features – We incorporated 15 geospatial habitat features, characterizing vegetation structure, plant species composition, and environmental variation, at varying resolutions (5-30 m), in our analysis of multiscale rattlesnake habitat use. Digital canopy height ("chm") and elevation models were developed from LiDAR data provided by the Ohio Geospatial Reference Program, collected in 2007 (; accessed 13 October 2014). The LiDAR data featured two returns pulse-1 with an average spacing and density of 1.27 m and 0.27 returns m-2, respectively. The canopy height and elevation models incorporated conventional methods, using bilinear interpolation, at 5 m resolution. From the elevation models, we developed layers on slope ("slope"), solar radiation ("sol_rad"), and the Beer’s transformed aspect ("beers"). In addition to canopy height represented in the "chm" column, we developed three penetration ratios (number of returns <2 m in height divided by the number of returns <50 m and <10 m, including returns <2m; and the number of returns <1 m divided by the number of returns <5 m, including returns <1 m) to characterize overstory ("ove"), midstory ("mid"), and understory ("und") vegetation density at 30 m resolution.

Woody plant composition was represented by three ordination axes within a gradient modeling approach developed from a separate study by B. Adams within the study area. Vegetation plot data, incorporating relative abundance profiles of trees and shrubs, including 99 woody plant taxa, were ordinated by non-metric multidimensional scaling (NMDS) and projected onto the landscape with terrain data and seasonal multispectral imagery provided by Landsat 8/OLI. Three subsequent floristic gradients synthesized moisture ("ndms1"), successional ("ndms2"), and elevational ("ndms3") variation among species responses within the vegetation data, at 30-m resolution. In addition, the plot data were used to develop geospatial layers on mean tree basal area ("rf_tba"; m2 ha-1) and density ("rf_tde"; stems ha-1).

From April 2017 to December 2017 we collected temperature data across the property using HOBO Pendant data loggers (Model # UA-001-08) set to one-hour logging intervals. We placed 150 loggers randomly across our study site by staking each logger in place at ground level on the north side of a tree to limit the influence of canopy structure on thermal data. We used these temperature data as a response variable in a linear mixed effects model with topographic and LiDAR-derived variables serving as predictive covariates. The final fitted model was used to predict the mean near-ground mid-summer temperature across our landscape ("meanT"). Finally, stand age ("stnd_age") was provided from the ODNR and a flowlines layer was transformed to a spatial grid (5 m) and used to summarize distances to the nearest stream ("strm_dist"; We extracted values for each covariate at all snake and random locations within different spatial scales of effect and centered and scaled all covariates prior to fitting statistical models.

Usage Notes

We modeled site use with generalized linear mixed effects Bernoulli models in a Bayesian framework using the brms package in R. The "dat_use_5m.csv" dataset is only used to run the general site use model (called in the "> Site Use Model" section of the R code). The four scale-specific datasets ("dat_5m.csv", "dat_25m.csv","dat_55m.csv","dat_105m.csv") are used to run the respective behavioral models.