Material and Methods
Study area
We studied snowshoe hares for three autumns (September
1st to December 1st of 2015, 2016,
and 2017) and four springs (March 1st to May
31st of 2015, 2016, 2017, and 2018) in southwestern
Yukon, Canada (Lat: 60.9 N, Long: -138.0 W). Snowshoe hares have been
monitored for over 40 years in this region (Krebs et al., 2018). Our
study area consists predominantly of white spruce (Picea glauca ),
trembling aspen (Populus tremuloides ) and balsam poplar
(Populus balsamifera ). Gray willow (Salix glauca ) and
dwarf birch (Betula glandulosa ) dominate the understory. The main
predators of snowshoe hares in this region include Canada lynx
(Lynx canadensis ), coyotes (Canis latrans ), goshawks
(Accipiter gentilis ), and great horned owls (Bubo
virginianus ) (Peers et al. 2020). Snowshoe hares went through
the increase, peak, and early decline phase of their population cycle
during our study period (Krebs et al. 2018).
Field methods
The study area was divided into three 35-ha snowshoe hare trapping
areas, located within ~ 8 km of each other (Peerset al. 2020). We captured snowshoe hares using Tomahawk
live-traps (Tomahawk Live Trap Co. Tomahawk, WI, USA) baited with
alfalfa and rabbit chow. Traps were set 30 minutes before sunset and
checked either three hours after sunset or at sunrise. We attached a
numbered ear tag to each hare to identify individuals on subsequent
recaptures, and we assessed coat-colour during each capture. To evaluate
coat colour, we examined hares from the front and sides and visually
estimated their percentage white coat to the nearest 5%. We later
binned coat colour in 10% white categories for analyses to account for
inter- and intra- observer ranking variability. We consider 10% bins as
reasonably precise given that intra- and inter-observer intraclass
correlation coefficients (ICC) for coat colour assessment were high
(ICC>0.9 in all cases, See Table S1 in Supplementary
Information). To monitor survival, we fit hares weighing >
1100g (n=347) with very high frequency (VHF) collars that were each
equipped with a mortality sensor (Model SOM2380, Wildlife Materials
Inc., USA, or Model MI-2M, Holohil, Canada, both < 27 ± 1 g).
We performed mortality checks of VHF collared hares almost daily, i.e.,
96.3% of checks occurred within 1 to 3 days. To monitor behaviour, we
also fit a subset (n =102) of VHF collared hares with an
accelerometer (model Axy3, 4 g, Technosmart, Rome, Italy).
Accelerometers measure force variation on three different axes and are
increasingly being used to infer behaviour in free-ranging animals
(Mikkelsen et al. 2019; Studd et al. 2019). Fully equipped
collars with both VHF and accelerometers had a total weight below 2.5%
of each individual’s body mass. Handling and collaring procedures were
approved by the University of Alberta Animal Care and Use Committee
(Protocol: AUP00001973).
We measured snow depth, snow cover, and temperature throughout our study
period. We measured snow depth on >60% of days at three
locations per trapping area, in relatively open forest, to the nearest
0.5 cm. Days with missing snow depth records were linearly interpolated
using the “zoo” function in the zoo package in R (Zeileis et
al. 2021). We measured snow cover by visually assessing daily landscape
photographs from three camera traps installed on each trapping area. We
calculated a combined average daily snow cover value to the nearest 10%
in our study region. We converted % snow cover to a binary type
variable above or below 60% snow cover (presence/absence) for the
autumn seasons, as there were very few instances when snow cover
estimates were between 0% and 100%. We measured temperature at least
six times a day on each trapping area using a minimum of 2 temperature
loggers (ibutton, DS1922L, Maxim Integrated, Whitewater, USA) to obtain
a single average daily temperature value for each trapping area.
Measuring coat colour mismatch
Coat colour mismatch was defined as the difference between hare percent
white (10% bins) and the daily percent snow cover (10% bins for both
autumn and spring). For all analyses, we treated mismatch as a binary
variable, defining mismatch as greater than 50% difference between hare
% white and snow cover (%). As such, mismatched hares were white
(> 50 % white) individuals in a snowless (< 50%
snow cover) environment. Considering that brown mismatched hares in a
snowy environment were rare (1% of trapping records), we did not
consider this type of mismatch in analyses. Although the threshold for
mismatch used in some previous studies is 60% contrast (Mills et al.
2013; Wilson et al. 2018), mismatch at this contrast threshold was rare
in our study region, i.e., in 11% of trapping records, so we used 50%
as our mismatch threshold to increase our sample size. That being said,
analyses using 40% or 60% thresholds for mismatch revealed similar
results (Tables S5, S6, S9, S10).
Effect of coat colour mismatch on survival
To evaluate the effect of coat colour mismatch on snowshoe hare
survival, we generated Cox’s proportional hazards (CPH) models (Cox &
Oakes 1984) with the “coxph” function in the survival package in R
(Therneau et al. 2021). The CPH model is a semi-parametric
approach used to analyze binary response data, in our case: alive or
dead (Sievert & Keith 1985). We monitored 347 hares and recorded 41
deaths over four springs and 34 deaths over three autumns. We excluded
mortality checks that exceeded seven days to limit the uncertainty in
the timing of death events (Murray & Bastille-Rousseau 2020). We
censored 15 individuals whose collars were removed before the end of the
study period and six individuals with permanently missing VHF signals.
We pooled data from different years, trapping areas, and sex, as
exploratory analysis indicated that none of those variables had a
significant effect on autumn or spring mortality risk (Table S2).
Considering that coat colour was assessed only during capture
opportunities (on average every 13.1 ± SD: 10.8 days per individual), we
assigned coat colour for each record in our survival analysis as the
nearest coat colour assessment completed in the field (average
difference of 4.95 ± SD: 3.70 days between telemetry check and
coat-colour assessment). We removed telemetry records where a coat
colour assessment within 14 days did not exist to ensure that coat
colour and derived mismatch values were an accurate representation of
each individual at the time of the telemetry check. Results from models
using survival records within 8 days of a coat-colour assessment were
qualitatively similar to those we obtained with our chosen 14-day
threshold (see Table S3).
We generated three competing CPH models for both autumn and spring. The
first model included snow cover and snow depth, based on prior evidence
of snow effects on hare survival (Meslow & Keith 1971; Peers et al.
2020). Our second model included those same snow variables in addition
to coat colour mismatch, our variable of interest. The third model was
the null (intercept-only) model. We used Akaike Information Criterion
for our model selection (Akaike 1974) and identified our top model based
on AICc (Burnham & Anderson 2002) with the package
AICcmodavg (Mazerolle 2019). We assessed multicollinearity in our top
model using the variance inflation factor (VIF) and ensured no variables
had VIF’s greater than 2. The proportionality assumption of CPH models,
which implies that the hazard ratio (HR; i.e., risk of death) is assumed
to be constant over time (Joshua Chen & Liu 2006), was met for our top
spring and autumn CPH model. Our results were not affected by
informative censoring, as we found qualitatively similar results for
both spring and autumn model coefficients when we treated censored
individuals as deaths (Murray & Bastille-Rousseau 2020) (Table S4).
Effect of coat colour mismatch on time spent foraging
To test our proposed mechanistic pathway, whereby white mismatched hares
experience reduced energetic requirements leading to reduced foraging
time (Balluffi-Fry et al. In Review; Sheriff et al.2009a), we used linear mixed-effects models using the “lmer” function
in the package lme4 (Bates et al. 2015). Daily time spent foraging
(minutes) was derived from tri-axial accelerometer data using
behavioural classifications previously developed in this hare population
(see Studd et al., 2019 for more information on classification methods).
Daily time spent foraging was classified over 4 second intervals at a
96% accuracy (Studd et al. 2019). We recorded 1505 daily
foraging records from 66 hares over the three autumns and 838 daily
foraging records from 44 hares over the four springs. Similar to our
survival analysis, we only kept foraging records that were within 14
days of a coat-colour assessment (average difference of 4.48 ± 3.51 (SD)
days). We reran our top foraging time models with data restricted to
daily foraging records that were within 8 days of a coat-colour
assessments instead to ensure that our results were not affected by this
14-day threshold, and obtained qualitatively similar results (Table S8).
To eliminate the potential of seasonal changes in foraging impacting our
results (Griffin et al. 2005), we restricted our data to only the
autumn and spring periods when snow cover was ≤ 50%, i.e., mismatch was
possible given our chosen threshold and therefore both matched and
mismatched individuals occurred simultaneously.
We generated four linear mixed-effects models per season to test for
differences in daily minutes spent foraging (our response for all
models) between matched brown hares and mismatched white hares and their
responses to changes in temperature. We included a random effect for
individual ID in all models to control for non-independence of data. We
included sex as a fixed factor in all spring models only, as exploratory
data analysis indicated that sex had a significant effect on time spent
foraging for spring but not autumn (Table S7). Furthermore, we included
year as a fixed effect in each model to account for potential effects of
yearly changes in predation risk on hare foraging behaviour (Shiratsuruet al. 2021). Our first model included two fixed effects,
temperature and year. Our second model included temperature, year, and
coat colour mismatch, and our third model included the same variables as
the second in addition to an interaction between mismatch and
temperature. Our fourth model was a null intercept-only model. We
checked model fit using marginal and conditional R- squared calculated
using the “r.squaredGLMM” function in the package MuMIn (Barton 2020),
according to Nakagawa et al. 2017. We used Akaike Information
Criterion (Akaike 1974) to rank our four competing models and identified
our top model in each season based on AICc (Burnham &
Anderson 2002).We completed all statistical analyses in R version 3.6.2
(2019) (R Core Team, 2019). We considered results where P ≤0.05 as
significant and reported all means with ± 1 standard error. The code
used for all analyses can be found online:
https://gitlab.com/joanie.kennah/chapter_1_kluane.