DOI: 10.5061/dryad.4b8gthtdm
title: Variation in behavior drives multiscale responses to habitat conditions in timber rattlesnakes (Crotalus horridus)

description: 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.
description: 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 &lt; 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 (http://ogrip.oit.ohio.gov/; accessed 13
 October 2014). The LiDAR data featured two returns
 pulse<sup>-1</sup> with an average spacing and density of 1.27
 m and 0.27 returns m<sup>-2</sup>, 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 &lt;2 m in height divided by the number of returns &lt;50
 m and &lt;10 m, including returns &lt;2m; and the number of
 returns &lt;1 m divided by the number of returns &lt;5 m,
 including returns &lt;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";
 m<sup>2</sup> ha<sup>-1</sup>) and density
 ("rf_tde"; stems ha<sup>-1</sup>).
 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"; https://nationalmap.gov/). 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.
description: 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 "&gt;
 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.