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.