Home range, site fidelity, and movements of timber rattlesnakes ( Crotalus horridus ) in west-central Illinois

Understanding the home range of imperiled reptiles is important to the design of conservation and recovery efforts. Despite numerous home range studies for the Threatened timber rattlesnake ( Crotalus horridus ), many have limited sample sizes or outdated analytical methods and only a single study has been undertaken in the central midwestern United States. We report on the home range size, site fidelity, and movements of C. horridus in west-central Illinois. Using VHF telemetry, we located 29 C. horridus (13 female, 16 male) over a 5-year period for a total of 51 annual records of the species’ locations and movements. We calculated annual home ranges for each snake per year using 99%,


Background
The home range of an animal was defined by Burt [1] as the geographic space "traversed by the individual in its normal activities of food gathering, mating, and caring for young".This definition has been well accepted since its conception and today is frequently characterized as the area encompassed by the animal's use of a geographic space (its "utilization distribution" or UD) [2].While home range estimation usually approximates an area surrounding observed or reported (e.g., through telemetry) locations, it also often involves some assumption (or estimation) of the pathway between these discrete locations.Methods used to connect these locations into pathways have evolved significantly since the home range was first defined [3,4].At its most simple, home range is defined by the use of the Minimum Convex Polygon (MCP) [5], where all the outer locations are connected to form a polygon and the home range is expressed as an area.Because MCP does not include information on the intensity of use within the polygon, it assumes equal use of the entire enclosed area, which can be an unrealistic appraisal of home range.However, due to the ease of calculating the MCP [2] and because it conservatively encompasses the entire range of an individual, it is still reported in the home range literature.Methods to calculate utilization distributions that include intensity-of-use within the range (e.g., Kernel Density Estimators) and probabilitystructured pathways based on knowledge of movement behavior (e.g., Brownian Bridge movement models) [6][7][8] have become more common and are usually considered a more accurate representation in visualizing space use or location choice by animals, despite increased computational complexity.
Evaluation of the home range of an animal can provide valuable insights into the ecological niche of a species, particularly when using utilization distributions.Utilization distributions delineate intensity-of-use areas that can be used to infer the habitat preference of individuals and provide a ready means to establish boundaries for microhabitat sampling efforts if it is assumed that increased presence is correlated with preference (e.g., [9][10][11][12][13][14][15][16][17].Delineation of high-use areas is particularly valuable when trying to understand critical areas for protection or enhancement with imperiled species [18][19][20][21].It can also be valuable in understanding habitat use (e.g., hunting or foraging areas) [22][23][24], how habitats might change annually or seasonally [21,[25][26][27][28] and in describing the site fidelity of individuals [29,30].
While the original definition of home range did not include a temporal component, one is usually provided [2], e.g., annual home range [31], seasonal home range [32], and internesting home range [28,33].Designating a temporal limit to the designation of home range provides the opportunity to compare animal habitat choice over time or under varying environmental conditions.
The imperiled timber rattlesnake (Crotalus horridus) is found throughout much of eastern North America, although their range has been reduced from its original distribution, and many local populations are extirpated or in decline [34,35].Crotalus horridus inhabits a wide diversity of habitats, including upland hardwood forests, mixed hardwood/pine forests, pine barrens and savannas, and coastal plain habitats, which include forest or forest/savanna mixes (Additional file 1: Table S1).Studies of C. horridus area use, home range size, and annual or daily movements are relatively common, particularly from the northeastern United States (Fig. 1; Additional file 1: Table S1), although many are unpublished graduate theses or suffer from limited sample sizes, low resolution of sampling, or outdated analytical methods which can affect home range estimates [36].Common findings among these studies include the observation that males have larger home ranges than females (e.g., [11,[36][37][38][39][40][41] and travel greater daily distances during their active period (e.g., summer) [36][37][38]42] and that gravid females have highly restricted home ranges that are much smaller than those of non-gravid females [11, 36-38, 40, 41].The larger male home range and daily movements are usually attributed to matesearching during mid-summer breeding periods [41].Despite the large number of studies on home range and movements of this species, there is little common agreement on home range sizes, movement distances and annual overlap in those high-use areas (see Additional file 1: Table S1).Some of this variance is likely due to limited sample sizes, low resolution of data collection [36] or outdated analytical methods and may also be the result of the wide range and habitats used by C. horridus.Conservation planning for this imperiled species is thus hampered, unless high-quality assessments are available within the range where such conservation actions are intended.
Only one study of the home range of C. horridus has been conducted in the central west part of the species range (Fig. 1), and no studies have been conducted in Illinois, where the species is listed as Threatened (https:// dnr.illin ois.gov/ conte nt/ dam/ soi/ en/ web/ dnr/ espb/ docum ents/ et-list-review-and-revis ion/ illin oisen dange redan dthre atene dspec ies.pdf ).Anderson [41] tracked 23 (12 males, 11 females) individual timber rattlesnakes for at least one active season in St. Louis County, Missouri and estimated the minimum and maximum "activity area" (which we interpret as a seasonal home range) for each snake using Minimum Convex Polygons.We calculated the mean home range for male, non-gravid and gravid female C. horridus using the data tables provided by Anderson: on average males had larger home ranges (x = 97.5 ha, SE = 24.8,n = 21) than non-gravid females (x = 12.1 ha, SE = 1.7, n = 22) and gravid females (x = 7.3, SE = 1.6, n = 11).Anderson [43] also reported that males moved greater daily distances during the mating season.While the Anderson studies [41,43] provide valuable insights into home range and area use of C. horridus in the west-central part of the species range, they were largely conducted in a highly mixed bottomland "glade" and upland forest mix, intermixed with developed/disturbed areas, and utilized only MCP as a measure of home range.Furthermore, like other studies their conclusions are limited by restricted sample size and tracking resolution.
In this study, we provide a second evaluation of home range and movement patterns for C. horridus in the west-central part of the species' range.In contrast to the Anderson studies [41,43], the environment in which our study was undertaken is composed primarily upland forest that has been undisturbed since the 1930s interspersed with a few old agricultural fields.Locations of active hibernacula have been known since the 1930s and protected since the 1950s.We also provide an assessment of site fidelity and annual overlap of high-use areas using updated analytical methods with Brownian Bridge utilization distributions and more refined data filtering techniques, and we include MCP for comparative purposes.S1 for further details on each study.Article citations marked with an asterisk denote unpublished reports or graduate theses.Geographic range of C. horridus is displayed for reference [60]

Study site and data collection
Our study was conducted at Principia College and surrounding properties in Jersey County, west-central Illinois (Fig. 1).The landscape is dominated by upland Oak-Hickory forest interspersed with occasional fields and human development.The site is bounded to the south by 75 m high limestone bluffs that run adjacent to the Mississippi River and are covered with a vegetational matrix of remnant hill prairies and woodlands.Crevices, scree slopes, and holes along the bluff front serve as hibernacula for C. horridus throughout the winter.Knowledge of these dens and their locations date back to the acquisition of the land for the College in the 1930s.
We captured C. horridus along ~ 3 km of bluff front during the species' fall ingress (September-October) and spring egress (April-May) periods or opportunistically throughout the active season (May-October).Captured snakes were weighed, measured (snout-vent length (SVL); tail-length (TL); subcaudal scale, and rattle segment counts) and sexed via cloacal probing, and a passive integrated transponder (PIT) tag (Avid Microchip ID Systems or BioMark PIT tag) was implanted subcutaneously in the distal one-third of the snakes' body.A subset of captured individuals was chosen for radiotelemetry based on size (i.e., ability to support a transmitter), sex (attempting to achieve an equal sex ratio), and reproductive status (gravid females identified by X-ray, were not implanted to prevent undue stress).Following the methods of Reinert and Cundall [44], temperature-sensitive radio transmitters (Holohil SI-2T, Carp, Ontario) were surgically implanted in the peritoneal cavities of snakes under inhalation anesthesia (isoflurane) by veterinarians at the Saint Louis Zoo.The transmitter weight was < 5% of a snake's body mass.Snakes were released at their original capture locations, usually within 1 week after surgery.We located C. horridus daily throughout their active season from their date of release, or spring egress in subsequent years, to fall ingress into hibernacula using a receiver (Advanced Telemetry Systems (ATS) model R410) and a three-element folding yagi antenna.At each location, we recorded the snake's geographic position using a handheld GPS device (Garmin eTrex 10 (2015-2019) or Bluetooth GPS (Garmin Glo) and Apple Ipad mini combination (2019), average accuracy of both systems ~ 3 m).

Home range estimation
All analyses were conducted using R Version 3.6.0[45].Maps were created in ArcMap 10.7 [46].We calculated annual home ranges for each snake per year using 99%, 95%, and 50% isopleths derived from Brownian Bridge utilization distributions (BBMM) [47].We also report 100% minimum convex polygons (MCP) [5] for comparison with other studies.BBMM requires two smoothing parameters, sig1 (related to the speed of an animal as calculated via its movement trajectory) and sig2 (related to the precision of GPS locations).Following the approach of Horne et al. [47], we calculated the maximum likelihood estimation of sig1 for each snake-year combination (using the 'liker' function in the R package 'adehabitat'), and we used the average GPS accuracy for all locations collected in the field (3 m) for sig2.We used area-asymptote curves to determine if enough fixes had been obtained for each snakeyear combination to estimate annual home range sizes [48,49].Each area-asymptote curve was estimated by removing all partial tracks (resulting from snake mortality or loss) from the data set and bootstrapping the MCP home range area of the remaining snake-year combinations by increments of 5 radiolocations (beginning at 5) for 500 iterations.We determined the point of asymptote by generating a line of best fit through the bootstrapped home ranges using nonlinear regression with a monomolecular growth function [50,51].We considered the number of fixes to be adequate if the home range size was > 95% of the asymptote.Snakes with non-asymptotic home ranges were removed from further analysis (but see methods on movement metrics below).
We examined the effects of sex, mass, and SVL on home range sizes independently for each of the four home range methods (99%, 95%, 50% BBMM isopleths, and 100% MCP) using linear mixed-effects models (LMM; R package 'nlme4') [52], with individual snake ID and year of tracking as crossed random intercepts in all candidate models.Although we had previously accounted for the effect of radio telemetry fixes by removing nonasymptotic home ranges in this analysis, we also included the number of fixes as a covariate in all models to control for varying sampling efforts.All candidate models were checked for normality of residuals using quantilequantile plots and Shapiro-Wilk tests, and homogeneity of residual variance using boxplots, Levene tests, and Pearson residual plots.While none of our models violated normality assumptions, all models which included sex exhibited heteroskedasticity; the home range size of males exhibited greater variance than females.Consequently, we refit all models which included sex with a within-group variance structure using the 'varIdent' function (R package 'nlme'), specifying sex as the grouping factor to allow for different variances across males and females.We then ranked candidate models, including a fully additive (global) and intercept-only (null) model, using Akaike's information criterion adjusted for small sample sizes (AIC c ; R package ' AICcmodavg' [53]).
We calculated the estimated marginal values for the top model in each candidate set (R package 'emmeans') and graphed (R package 'ggplot2' [54]) the subsequent values and 95% confidence intervals.Parameters with 95% CIs not overlapping zero were considered significant influences of annual home range size.

Movement metrics
We calculated commonly reported movement metrics, including maximum displacement (i.e., distance) from hibernacula and mean daily movement distance.Daily displacement distances were calculated for each snakeyear combination as the Euclidian (straight-line) distance between a snake's recorded location and its hibernacula for each day of its tracking duration.Daily movement distances for each snake-year combination were calculated as the Euclidean distance between a snake's recorded location on each day and the previous day's location.If multiple days had passed between successive locations, the distance moved was divided by the number of days between locations.Like the home range analysis, we examined the effects of sex, mass, and SVL on maximum displacement distance from hibernacula and average daily movement distances using linear mixed-effects models with the same model specification and AIC c structure as before.All models satisfied assumption tests (see above).The top-ranked AIC c model(s) were then graphed to display both predicted values and 95% CIs.
We also examined temporal changes in displacement distances and daily distances traveled throughout the active season by plotting their daily arithmetic averages across the tracking period against the day of year.To avoid the problem of pseudo-replication, we employed a resampling approach to calculate daily means and standard errors of both movement metrics.Specifically, for each day of the year, we randomly sampled one value per snake without replacement and calculated the average of the resulting samples.We repeated this process for 5000 iterations and then calculated the average and standard error of the resulting mean distribution.Bootstrapping ensured that each snake was only represented once in any given calculation while also allowing snakes with multiple years of the data to be represented.Unlike our home range analyses in which snake-year combinations that failed to reach asymptote were removed, all snakes were included in the daily movement analysis.These daily averages (± standard error) were displayed as 7-day moving averages grouped by sex to allow for examination of general phenological trends.

Home range overlap and site fidelity
We calculated intra-individual home range overlap (i.e., multi-year site fidelity) for C. horridus with multiple annual home ranges that met the previously described asymptotic criteria using several methods.Following the recommendations of Fieberg and Kochanny [55], we estimated the overlap for the 99%, 95%, and 50% BBMMderived home ranges using the utilization distribution overlap index (UDOI) and Bhattacharyya's affinity (BA: [56]) both of which incorporate the use-intensity aspect of the utilization distribution in overlap estimations.We also calculated simple isopleth overlaps for the BBMM isopleths and 100% MCP estimates by dividing the area of intersection between two given home range isopleths by the total combined area of both isopleths.For each overlap metric, we calculated the pairwise overlap between every intra-individual snake-year home range combination.To account for pseudo-replication (multiple overlaps per snake), we calculated the estimated mean values and 95% CIs of each overlap method by fitting Beta Generalized Linear Mixed Models (r package 'glmmTMB' [57]) with Snake ID as a random effect.We specified a sex by isopleth size interaction as the fixed effects for models examining BBMM home range overlap (99%, 95%, and 50%) to account for differing overlap values across isopleth size.We specified only sex for the model examining MCP home range overlap.We graphed the resulting model intercept estimates representing the mean of each overlap metric across sex.

Results
We radio-tracked 29 individual C. horridus (13 female, 16 male) for one or more years between 2015 and 2019, accumulating a total of 51 annual records (Table 1).Nine of the 51 annual records were partial tracks (due to snake mortality or signal loss), and another six failed to reach the home range asymptote, suggesting that they had an inadequate number of radio-location fixes to estimate annual home range sizes and were subsequently removed from further analysis (Table 1 ).We note that females identified here as gravid were not identified as gravid at radio tag implantation, but likely became gravid after that implantation or embryo development was not apparent by X-ray.

Home range estimation
Prior to analyzing factors that might account for variation in home range size, we removed gravid female home ranges (n = 2) from the data set due to small sample sizes and their uniquely distinct home range size and characteristics.The candidate model that included the independent variables sex and number of fixes was the top AICc-ranked model for all home range methods (Additional file 1: Table S3).As indicated via 95% confidence intervals, sex was a significant influencer of home range size for all methods of home range estimates (Additional file 1: Table S4), with male C. horridus exhibiting significantly larger home ranges than non-gravid females (Fig. 2A-D).Conversely, the number of fixes was insignificant in all top models.The predicted (marginal) home range sizes for male and non-gravid C. horridus when controlling for the number of fixes were 88.72 Ha (CI 63. 41

Movement metrics
Similar to the home range analysis, the candidate model including the independent variables sex and number of fixes was the top AICc-ranked model for both daily distance traveled and maximum displacement distance from hibernacula (Additional file 1: Table S3).Sex was a significant influencer in both movement metrics (Additional file 1: Table S4), whereas the number of fixes was insignificant.The estimated daily distance traveled was significantly greater for males (mean = 57.25 m/day, CI 49.06-65.43)than females (mean = 27.55 m/day, CI 18.99-36.12)(Fig. 2E), with a marked increase in male movements from late June (~ day of year 176) to mid-August (~ day of year 231) (Fig. 4B).Similarly, maximum Fig. 2 Marginal effects of sex (F: female; M: male) on the home range area and movement estimates (controlling for the number of GPS fixes) for timber rattlesnakes (Crotalus horridus) as predicted by top AIC c -ranked linear mixed-effects models.Subplots represent predicted home range areas estimated using Brownian Bridge Movement Models (BBMM) at 99% (A), 95% (B), and 50% (C) isopleths, and minimum convex polygon (MCP; D), and predicted estimates of mean daily distance travelled (E), and maximum displacement distance from over-wintering hibernacula (F).Plots indicate predicted values (solid black dots) ± 95% confidence intervals (CIs).Data derived from 36 snake-year combinations from 21 individuals' radio-tracked from 2015 to 2019 in Jersey County, Illinois.See Additional file 1: Table S3 and 4 for a AIC c table and parameter estimates and on average, males were located further from their hibernacula throughout the entirety of their active season (Fig. 4A).

Discussion
We used VHF telemetry to record the daily locations of 29 free-ranging C. horridus during their summer activity period from 2014-2019 for a total of 51 activity periods at a site in west-central Illinois.After filtering the dataset using area-asymptote curves to remove those with inadequate sample sizes, 36 of the tracking records were sufficient for home range estimation.The home ranges for males were larger than those for non-gravid females, and both were larger than the home range of gravid females.Due to the wide range of analytical methods used by studies of the home range for this species, a direct comparison of our home range size to others is difficult and standardization of study methods across all such studies is needed.Nonetheless, we note that the home range of our males and non-gravid females fall close to the mean of all home range utilization distributions (85.5 Ha and 36.2Ha) reported for the species in the published and unpublished scientific literature (Additional file 1: Table S1).Dispersal distances from hibernacula and daily travel distances were also greater for males (Figs. 2 and  4).Such results are consistent with other studies (e.g., [11,[36][37][38][39][40][41]43]), where it is proposed that the search behavior of males looking for females during the mating season results in larger home ranges.Males in our study showed a distinct increase in movements between ~ day 185 (4 July) and ~ day 235 (22 August) each summer (Fig. 4), which encompasses the mating season for this location [58].Even though sex was an important variable affecting home range size, there was still a high amount of unexplained variance, particularly for males (Fig. 2).While we did not have enough statistical power to explore interactions between sex and other variables, we suspect there are significant interactions between sex and SVL or sex and mass [36].Males are the active participants in mate-searching, and thus larger males (in mass or SVL) will probably exhibit larger home ranges as they search for females, whereas non-searching immature males will exhibit smaller home ranges.This likely explains the larger variation in home range size exhibited by males in our study.Conversely, the home range size of females was less variable, probably because they are inactive participants in mate-searching and thus range size varies less.Sex and size are likely correlated with the sexual size dimorphism of this species with males being typically larger than females, so careful analysis of such interaction is warranted.
When evaluating general trends in dispersal from hibernacula, two patterns emerged (Fig. 3).Females move directly away from their hibernacula to their high-use areas, and their home ranges appear very linear.While males also move directly to high-use areas, their use areas are much broader, likely due to mate-searching in mid-summer, so their overall home ranges look less linear than those of female rattlesnakes.Second, all snakes from specific hibernacula tend to move in a similar dispersal direction (Fig. 3).The snakes from hibernacula 1 (located ~ 0.74 km east of hibernacula 2) dispersed from their hibernacula toward the northwest, and snakes from hibernacula 2 moved directly north or slightly NE.The difference in dispersal direction and lack of home range overlap between snakes dispersing from different hibernacula might be due to potential barriers to dispersal, which in this case are old agricultural fields located north An alternate and possibly complimentary explanation for the dispersal pattern following fixed directions from their hibernacula might be that dispersal direction and the subsequent high-use area are linked to timber rattlesnake tendency for natal homing.That offspring return to their natal hibernacula has been theorized as an explanation for genetic differentiation between hibernacula [41,[58][59][60].If timber rattlesnakes have strong natal homing and follow their mothers into high-use areas on first departure from their birth areas or during their first spring, it is conceivable that they may imprint on those Fig. 5 Average annual home range overlap indices (± CI) for 11 timber rattlesnakes (Crotalus horridus; n = 20 home range overlaps) tracked for two or more years in Jersey County Il.Overlap of home ranges calculated using the Brownian Bridge Movement Model (BBMM; for 99%, 95% and 50% isopleths) were quantified using Bhattacharyya's affinity (BA; A), the Utilization Distribution Overlap Index (UDOI; B), and simple isopleth overlaps (i.e., dividing the area of intersection between two given home range isopleths by the total combined area of both isopleths; C).The latter metric was also used to calculate overlap of home ranges calculated via Minimum Convex Polygons (MCP; D) high-use areas, much as they seem to imprint on their natal hibernacula.That such areas have been successful for their mothers (as evidenced by successful reproduction), might imply potential success for the offspring and thus promote this behavior within the species.Evaluating the maternal origins of foraging snakes (e.g., using MtDNA) in distinct foraging areas could yield valuable insights into how this species chooses foraging locations.
Irrespective of the reason for dispersal direction, individual snakes tended to return to the same use areas each year.The average home range overlap at 99% BBMM was 0.49 and 0.78 using UDOI overlap, and both values imply relatively high overlap between years, particularly when considering all the areas a snake may choose to inhabit during the active season.Similar results in high-use area overlap have been shown by Nordberg et al. [61] in Tennessee.These significant levels of return to the same high-use areas may partially explain why relocation of timber rattlesnakes show reduced fitness [38].Once established on a specific use area, timber rattlesnakes may lack adaptability to forage in alternate areas and relocated snakes do not thrive.Nonetheless it is apparent from the reduced overlap in the 50% UD that there is some plasticity in choice of high-use areas which may reflect flexibility in the hunting strategies of this apex predator.

Conclusions
Our study of the home range and active-season movements of timber rattlesnakes provide supporting information on the sexually dimorphic home range size and movement patterns of the species in the central part of the species range.We believe that the high quality of our estimations of home range and use area overlap can provide useful parameters for the implementation of conservation actions in this portion of this imperiled species range.Our observations of high levels of annual home range overlap and possible segregation of high-use areas by hibernacula also provide further insights into the biology of this important forest predator as well as supporting the idea that relocation of timber rattlesnakes is not a viable conservation protocol.

Fig. 1
Fig.1Location of the current study (yellow circle) and other studies which report timber rattlesnake (Crotalus horridus) home range statistics (red circles).See Additional file 1: TableS1for further details on each study.Article citations marked with an asterisk denote unpublished reports or graduate theses.Geographic range of C. horridus is displayed for reference[60]

Fig. 3
Fig. 3 Home range maps for 36 snake-year combinations from 21 timber rattlesnakes (Crotalus horridus) tracked between 2015 and 2019 in Jersey County, Il.Each map displays: GPS locations (fixes) from which home range estimates were calculated (black dots); 99% (red), 95% (yellow), and 50% (orange) BBMM home range estimates; MCP home range estimates (light grey); and location of hibernacula (blue star).Annotations indicate hibernacula ID (bottom right); snake ID and year of track (top left); and sex and reproductive condition (M: male, F: female, G: gravid; top right).Maps ordered by hibernacula, sex, and year

Table 1
Summary statistics for 51 timber rattlesnake (Crotalus horridus) snake-year combinations from 31 individuals (13 female, 16 male) tracked between 2015 and 2019 at Principia College and surrounding properties in Jersey County, Illinois Columns represent: snake ID (ID); year of track (Year); sex and reproductive condition (sex; M: male, F: female), G: gravid); hibernacula (Den ID); snake mass (Mass); snout-vent length (SVL); start (Track start) and end (Track end) date of track; track duration (Track length); number of GPS coordinates (# fixes); average number of days between fixes (Days / fixes); minimum number of fixes required for home range asymptote (Min.fixes); proportion of home range (Prop.HR); mean daily distance (Daily dist.);maximum distance/displacement from hibernacula (Max.disp.);total planimetric distance moved between fixes (Total dist.

Table 2
Grand means (± SE) of snake mass (Mass), snout-vent length (SVL), radio-tracking duration (Track length), number of GPS locations (# fixes), days between fixes (Days/fixes), and minimum fixes required for home range asymptote (Min.fixes) for 36 snakeyear combinations from 21 individual timber rattlesnakes (Crotalus horridus), tracked for one or more years in Jersey County Illinois, during the active seasons of 2015-2019 Data summaries include multiple entries for individuals tracked for one or more years