https://cpw.cvlcollections.org/files/original/4aeb43ecac4e92751417b9a77893956c.pdf 5e07d0694764da3e41c69d0c78071ac1 PDF Text Text The research in this publication was partially or fully funded by Colorado Parks and Wildlife. Dan Prenzlow, Director, Colorado Parks and Wildlife • Parks and Wildlife Commission: Marvin McDaniel, Chair • Carrie Besnette Hauser, Vice-Chair Marie Haskett, Secretary • Taishya Adams • Betsy Blecha • Charles Garcia • Dallas May • Duke Phillips, IV • Luke B. Schafer • James Jay Tutchton • Eden Vardy �Received: 22 June 2019 | Revised: 23 October 2019 | Accepted: 20 November 2019 DOI: 10.1111/geb.13049 RE SE ARCH PAPER Local climate determines vulnerability to camouflage mismatch in snowshoe hares Marketa Zimova1 | Alexej P. K. Sirén2,3 | Joshua J. Nowak1 | 3,4 5 2,3 Alexander M. Bryan | Jacob S. Ivan | Toni Lyn Morelli | 1 6 1,7 Skyler L. Suhrer | Jesse Whittington | L. Scott Mills 1 Wildlife Biology Program, University of Montana, Missoula, Montana Abstract 2 Aim: Phenological mismatches, when life-events become mistimed with optimal Department of Environmental Conservation, University of Massachusetts, Amherst, Massachusetts 3 U.S. Geological Survey, Northeast Climate Adaptation Science Center, Amherst, Massachusetts 4 Integrated Environmental Data, LLC, Berne, New York 5 Colorado Parks and Wildlife, Fort Collins, Colorado 6 Parks Canada, Banff National Park Resource Conservation, Banff, Alberta, Canada 7 Wildlife Biology Program and Office of Research and Creative Scholarship, University of Montana, Missoula, Montana Correspondence Marketa Zimova, Wildlife Biology Program, University of Montana, Missoula, MT 59812, USA. Email: mzimovaa@umich.edu Present address Marketa Zimova, School for Environment and Sustainability, University of Michigan, Ann Arbor, Michigan Funding information University of Montana; U.S. Geological Survey, Grant/Award Number: Department of the Interior Northeast Climate Adaptation Science Center and Department of the Interior Southeast Climate Adaptation Science Center; Parks Canada; North Carolina State University; Vermont Fish and Wildlife Department, Grant/Award Number: E-1-25 grant for Investigations and Population Rec; Division of Environmental Biology, Grant/Award Number: National Science Foundation Division of Environmen and National Science Foundation EPSCoR Award No. 17362; New Hampshire Fish and Game Department, Grant/Award Global Ecol Biogeogr. 2020;29:503–515. environmental conditions, have become increasingly common under climate change. Population-level susceptibility to mismatches depends on how phenology and phenotypic plasticity vary across a species’ distributional range. Here, we quantify the environmental drivers of colour moult phenology, phenotypic plasticity, and the extent of phenological mismatch in seasonal camouflage to assess vulnerability to mismatch in a common North American mammal. Location: North America. Time period: 2010–2017. Major taxa studied: Snowshoe hare (Lepus americanus). Methods: We used > 5,500 by-catch photographs of snowshoe hares from 448 remote camera trap sites at three independent study areas. To quantify moult phenology and phenotypic plasticity, we used multinomial logistic regression models that incorporated geospatial and high-resolution climate data. We estimated occurrence of camouflage mismatch between hares’ coat colour and the presence and absence of snow over 7 years of monitoring. Results: Spatial and temporal variation in moult phenology depended on local climate conditions more so than on latitude. First, hares in colder, snowier areas moulted earlier in the fall and later in the spring. Next, hares exhibited phenotypic plasticity in moult phenology in response to annual variation in temperature and snow duration, especially in the spring. Finally, the occurrence of camouflage mismatch varied in space and time; white hares on dark, snowless background occurred primarily during low-snow years in regions characterized by shallow, short-lasting snowpack. Main conclusions: Long-term climate and annual variation in snow and temperature determine coat colour moult phenology in snowshoe hares. In most areas, climate change leads to shorter snow seasons, but the occurrence of camouflage mismatch varies across the species’ range. Our results underscore the population-specific susceptibility to climate change-induced stressors and the necessity to understand this variation to prioritize the populations most vulnerable under global environmental change. wileyonlinelibrary.com/journal/geb�© 2019 John Wiley & Sons Ltd | 503 �504 | Number: T-2-3R grant for Nongame Species Monitoring and Ma; National Science Foundation, Grant/Award Number: 0841884 and 1907022; National Science Foundation, Grant/Award Number: EPCSoR 1736249 (OIA-1736249); Colorado Species Conservation Trust Fund ZIMOVA et al. KEYWORDS adaptation, camouflage mismatch, climate change, latitudinal gradient, phenological mismatch, phenotypic plasticity, range edge, snow, snowshoe hares Editor: Patricia Morellato 1 | I NTRO D U C TI O N consequences for susceptibility to phenological mismatch on both the local population and broader species levels. As a result of anthropogenic climate change, plant and animal popu- Seasonal coat colour moult, a key phenological trait, has received lations are increasingly confronting environmental conditions differ- increased attention as a trait shaped by climate. Across the Northern ent from the ones to which they are adapted. Organisms occupying Hemisphere, 21 species of mammals and birds change coat or plum- seasonal environments have evolved mechanisms for the timing of age colour from brown in the summer to white in the winter to match their life cycles (i.e., phenology) to match with optimal environmen- snow-covered landscapes (Mills et al., 2018; Zimova et al., 2018). As tal conditions and resources at their location (Bradshaw & Holzapfel, with other phenological traits, photoperiod serves as the principal cue 2007; Williams, Henry, & Sinclair, 2015). When phenology and fa- for moult phenology, with some evidence that year-to-year variation vourable environmental conditions become asynchronized, organ- in winter weather modulates the progression of moult (Hofman, 2004; isms can suffer negative fitness costs (Both, Bouwhuis, Lessells, & Zimova et al., 2018). However, decreasing duration of snow cover due Visser, 2006; Lane, Kruuk, Charmantier, Murie, & Dobson, 2012; to climate change (Choi, Robinson, & Kang, 2010; Kunkel et al., 2016; Post & Forchhammer, 2008; Senner, Stager, & Sandercock, 2017; Vaughan et al., 2013) may result in phenological mismatch, whereby Zimova, Mills, & Nowak, 2016). Such phenological mismatches winter white animals become colour mismatched against dark, snow- are becoming increasingly common under climate change (Cohen, less backgrounds (Mills et al., 2013). Field studies indicate that mis- Lajeunesse, & Rohr, 2018; Parmesan & Yohe, 2003; Thackeray et al., match in seasonal coat colour and snow presence or absence leads to 2010), which in the absence of adaptive responses could lead to pop- high fitness costs due to increased predator-induced mortality (Atmeh, ulation declines and local extinctions (Visser & Gienapp, 2019). Thus, Andruszkiewicz, & Zub, 2018; Zimova et al., 2016) and may have already there is a pressing need to understand the degree to which climate contributed to range contractions for several species including Lagopus change leads to phenological mismatches and the capacity of wild and Lepus spp. (Diefenbach, Rathbun, Vreeland, Grove, & Kanapaux, populations to withstand attendant fitness costs. 2016; Imperio, Bionda, Viterbi, & Provenzale, 2013; Pedersen, Odden, Understanding the variation in phenology and its environmen- & Pedersen, 2017; Sultaire et al., 2016). tal drivers is fundamental for assessing current and future species’ The snowshoe hare (Lepus americanus), a key prey species of the vulnerability to phenological mismatches. For the majority of traits boreal forest of North America (Krebs, Boonstra, Boutin, & Sinclair, in plants and animals in temperate regions, photoperiod serves as 2001), exhibits seasonal colour moults in the majority of its range the principal cue for phenology, with temperature and other envi- (Gigliotti, Diefenbach, & Sheriff, 2017; Mills et al., 2018; Nagorsen, ronmental factors exerting lesser influence (Bradshaw & Holzapfel, 1983). Because of their broad distribution (Figure 1a), hares inhabit a 2007; Hofman, 2004). Because photoperiod, latitude and climate co- large range of environmental conditions, making them an ideal spe- vary across most species’ ranges, variation in phenology is often dis- cies for investigating variation in moult phenology and camouflage tributed along latitudinal gradients. Northern populations typically mismatch. In the only two areas where moults have been recently in- experience harsher and colder climates that correspond with later vestigated in relation to climate change (i.e., Montana and Wisconsin, initiation of spring and earlier initiation of fall events compared to USA), phenotypic plasticity is not sufficient to prevent camouflage southern populations (Bradshaw & Holzapfel, 2007; Hut, Paolucci, mismatch (e.g., Mills et al., 2013; Wilson, Shipley, Zuckerberg, Peery, Dor, Kyriacou, & Daan, 2013). However, latitude may not be a reli- & Pauli, 2018). For example, hares in Montana experience about a able predictor of phenology for two reasons. First, the covariance week of mismatch annually whereby hares are in the wrong coat co- between latitude and climate is imperfect when other geographic lour in relation to their background (i.e., white hares on snowless factors such as elevation also affect climatic conditions (Chuine, background or brown hares on snow; Mills et al., 2013; Zimova, 2010; Visser, Caro, Oers, Schaper, & Helm, 2010). Secondly, species Mills, Lukacs, & Mitchell, 2014; Zimova et al., 2016). Because hares show year-to-year in situ variation in phenological traits in response rely heavily on their camouflage for survival, mismatch has strong to annual variation in climate. This temporal phenological variation negative fitness costs (i.e., 7–12% reduced weekly survival) that in is referred to as ‘population-level’ phenotypic plasticity, and is dif- the absence of evolutionary shifts and under projected snow de- ferent from between-individual level phenotypic plasticity (Gienapp clines would be sufficient to cause steep population declines and & Brommer, 2014; Phillimore, Hadfield, Jones, & Smithers, 2010). local extinction (Wilson et al., 2018; Zimova et al., 2016). To date, Overall, both spatial and temporal variation in phenology may have no study has evaluated the phenological drivers and camouflage �| ZIMOVA et al. F I G U R E 1 Camera site locations and snowshoe hare moult phenologies and moult dates in the Canada, Colorado and New England study areas. (a) Snowshoe hare range (downloaded from www.iucnr​ edlist.org) is coloured and shaded by the mean annual number of snow days (Dietz, Kuenzer, & Dech, 2015). Coloured points represent the 448 remote-camera trap sites. (b) Bold lines depict predicted probabilities of being in the final colour category (white in the fall, brown in the spring) over time. The dashed lines show 95% credible intervals. The horizontal dashed lines at .90 intersect with population means to identify moult completion dates. Population mean moult initiation and completion dates are depicted as a date range in the bottom right corners, with the completion dates in bold. The predicted probabilities and dates were estimated for each season and population based on the model without covariates [Colour figure can be viewed at wileyonlinelibrary.com] 505 (a) (b) mismatch across the heterogeneous snowshoe hare range in a uni- encompass wide latitudinal, altitudinal and habitat variation across fied analytical framework. the distributional range of snowshoe hares (Table 1). Our northern- To understand snowshoe hare vulnerability to camouflage mis- most study area in the Canadian Rockies included three national match, we quantified the spatial and temporal variation in colour moult parks (Banff, Yoho and Kootenay National Parks) and was charac- phenology, phenotypic plasticity, and the occurrence of camouflage terized by rugged, forested mountainous regions with a long snow mismatch across three disjunct, climatically and geographically distinct season. This area is located within a homogeneous boreal habitat at populations. First, we hypothesized that the spatial variation in moult the core of the snowshoe hare distribution range (Cheng, Hodges, phenology would be determined by latitude and local climate. We Melo-Ferreira, Alves, & Mills, 2014). The southernmost study area tested two predictions: (a) populations in more northern sites moult was located in the San Juan Mountains of south-west Colorado; an earlier in the fall and later in the spring, (b) populations in colder and isolated patch of high-elevation southern boreal forest near hares’ snowier sites moult earlier in the fall and later in the spring. Second, southern range boundary. The northern New England study area, we hypothesized that hares exhibit temperature- and snow-mediated also near the southern edge of hare distribution, encompassed por- phenotypic plasticity in moult phenology and we predicted that moults tions of the Green Mountains in Vermont and the White Mountains occur earlier in fall and later in spring during colder and or snowier in New Hampshire. This area stretched across a transition zone be- years. Third, we quantified the occurrence of camouflage mismatch for tween the northern hardwood and boreal forests and had, on av- each population in spring and fall and assessed which snowshoe hare erage, the mildest climate of the three areas. Major hare predators populations may be the most vulnerable to camouflage mismatch now in all areas included coyotes (Canis latrans), red fox (Vulpes vulpes), and under future climate change. Canada lynx (Lynx canadensis), bobcats (Lynx rufus), martens (Martes spp.), weasels (Mustela spp.), northern goshawks (Accipter gentilis) and 2 | M ATE R I A L S A N D M E TH O DS 2.1 | Study areas great-horned owls (Bubo virginianus). The study area in Canada also had wolves (Canis lupus) and grizzly bears (Ursus arctos). 2.2 | Camera trapping design Our analysis integrated three regional monitoring studies in North America: the Canadian Rockies, the San Juan Mountains in Colorado, We obtained our snowshoe hare phenology data as by-catch pho- and northern New England (Figure 1a). Together, these areas tos from three long-term independent studies of forest carnivores. �| 125.01 (12.49) 155.70 (11.07) 168.69 (18.18) Snow season (days) ZIMOVA et al. A combination of motion-triggered camera models was used across the sites, but all produced comparable high-quality images during the day and night. The camera spacing differed between regions, 8.35 (1.84) 8.03 (1.71) 5.17 (2.33) tmax (°C) but at any given time the minimum distance between cameras was > 1 km; grid sizes varied by study area (10 km × 10 km in Canada, 4 km × 4 km in Colorado and 2 km × 2 km in New England). Each camera site was stationary and remained within a grid cell for a minimum −3.44 (0.93) −7.24 (1.28) −6.62 (1.68) tmin (°C) duration of 1 year. We combined all observations for the cameras that were moved between years within < 1 km distance and assigned them an average latitude, longitude and elevation. The 1-km threshold was chosen because it far exceeds average individual hare 1,122 1,705 921 105 52.07 (7.30) 10.61 (1.63) Coat colour was visually estimated from the by-catch hare photographs by one observer, following a standardized protocol for colour classification and image quality filtering (Zimova et al., under review). Briefly, hares were classified as (a) white when > 90% of the body (excluding belly and feet) was white, (b) brown when < 10% of the body was white, and (c) moulting for all other instances. As per previously developed classification methods (Zimova et al., 2014), we excluded the feet and belly colour as these stay white all year in all three areas (Grange, 1932). Because hares cannot be individually identified from photos, we considered 0.25 (0.89) 183 72.83 (6.91) 9.70 (1.78) −4.68 (1.09) 91 93.08 (10.93) 5.22 (1.86) tmax (°C) tmin (°C) 1999; Ivan, 2011; Mills et al., 2005). 2.3 | Coat colour monitoring −5.41 (1.33) Camera sites (n) Snow season (days) Spring Obs. (n) seasonal movement and home range size at all three areas (Hodges, photos spaced by > 24 hr as independent events, unless different individuals were distinguishable (e.g., two hares in one photo, a hare in a different moult stage recorded within 24 hr). 486 322 November 2016, yielding 1,888 independent coat colour observa967 Obs. (n) In Canada, 121 cameras operated from September 2011 to tions. In Colorado, 206 cameras operated from September 2010 to 65 from January 2014 to January 2018, yielding 1,608 independent observations. 627.18 (284.59) 110 3,200.79 (215.46) 98 independent observations. In New England, 121 cameras operated 1,850.98 (266.20) Elevation (m) Camera sites (n) June 2011 and from January 2014 to August 2017, yielding 2,027 2.4 | Statistical analysis We used R (version 3.4.3; R Core Team, 2016) for all statistical 44.54 (0.51) 37.63 (0.30) 51.38 (0.33) Latitude (degrees) analyses. 2.4.1 | Climate variables To characterize annual and long-term climate conditions at each New England Colorado Canada camera site, we prepared a set of temperature and snow cover variStudy area Fall TA B L E 1 Geospatial and long-term climate details regarding the camera trap networks in the study areas. Number of camera sites and independent coat colour observations are given separately for each season and area (across all years and areas the number of independent camera sites totaled 448). Other metrics include mean values across all camera sites within a study area with standard deviation in parentheses. Long-term mean minimum (tmin) and maximum (tmax) temperature and snow season duration are based on 1980–2009 period 506 ables relevant to moult phenology (Mills et al., 2018; Zimova et al., 2018). Annual minimum (tmin) and maximum (tmax) temperature during spring and fall were calculated for each year of monitoring (2010–2017) at each camera site using daily gridded 1 km × 1 km �| ZIMOVA et al. 507 resolution data (Daymet; Thornton et al., 2018). We used the same hares initiated and completed their moults as ‘initiation’ and ‘com- dataset to calculate mean seasonal tmin and tmax during the 30- pletion’ dates at each area. Fall start was specified as the first Julian year period 1980–2009 to describe long-term climate at each cam- day when mean pbrown < .9 and end date when mean pwhite > .9; the era site. The seasons were defined as spring (1 March–31 May) and opposite condition was used to estimate the spring dates (i.e., start fall (1 September–30 November), because at all areas the majority of when pwhite < .9 and end when pbrown > .9). moulting occurs during those months. Next, to test the effect of local environmental covariates on phe- We used modelled snow water equivalent (SWE) to quantify the nology, we combined colour observations from all years and popula- duration of the continuous snow season and the total number of snow tions in one dataset and constructed a set of univariate models. Each days for the spring and fall. For the years of monitoring (2010–2017), model included a single fixed effect of an environmental covariate we used daily gridded 1 km × 1 km resolution data by Snow Data β2i on the probability of being brown, white or moulting. The envi- Assimilation (hereafter SNODAS: Barrett, 2003). Because SNODAS ronmental covariates were latitude, elevation, and the 30-year long- is unavailable prior to 2000, we used daily gridded 6 km × 6 km reso- term temperature and snow conditions (i.e., tmin and tmax, duration lution data (Livneh et al., 2015) to describe the mean long-term snow of snow season) in spring and fall during each season at each cam- conditions during the 30-year period (1980–2009). The duration of era site. We used univariate models to avoid problems associated spring snow season each year was calculated as the number of days with multicollinearity as most environmental covariates were highly between 1 January and snowmelt date, that is, the first day when correlated (Pearson correlation coefficients > |.60|; Supporting snow is absent (SWE = 0) at a camera site for a minimum of 7 days. Information Table S1). To facilitate comparisons between models, all Fall snow season was estimated as the number of days following the covariates were standardized to a mean of 0 and SD of 1. summer (i.e., the longest snowless period at each camera site each The resulting β coefficients represented the increase in the prob- year), between snow onset date (i.e., first day when SWE > 0 for a ability of being in a certain colour category on the multinomial-logit minimum of 7 days) and 31 December. We used the 7-day buffers scale for every one-unit (SD) change in the covariate. We considered to discard spurious early spring snowmelts followed by further ex- a covariate to have a significant effect on moult phenology when the tended snow season and to account for spurious snow flurries in the resulting β coefficients’ 95% credible intervals (CRIs) did not over- fall. The total number of snow days was calculated as the sum of days lap zero. Because we were not interested in quantification of the with SWE > 0 from 1 January to 31 August for the spring and from 1 effect size per se, but rather on the direction and the relative sizes September to 31 December for the fall. of the different covariates, we did not convert the coefficients to normal scale. We primarily focused on the covariate effects on probability of the season’s final colour (i.e., β2white in the fall, β2brown in 2.5 | Moult phenology the spring) not the initial colour or moulting colour category, in part to simplify the reporting of results. Furthermore, because photope- We used a hierarchical multinomial logistic regression analysis riod is known to strongly control the hormonal cascade that triggers within a Bayesian framework to describe moult phenology and its the moult, we expected the effects of climate to be more appar- phenotypic plasticity. For all models, we estimated the probability of ent in the final rather than the initial stage of the moult (i.e., follicle a hare being in colour category i at a camera site j on a Julian day d as: stimulation and hair growth initiation versus the appearance of the newly grown hairs and shedding of the old hairs; Zimova et al., 2018). eαi +β1i ∗d+si,j p (y = i) = ∑i−1 αi +β1i ∗d+si,j 1 + k=1 e (1) Finally, in most cases, the significance or absolute effect size of initial and final colour probabilities were similar (all β coefficients shown in Supporting Information Tables S2 and S3). Coat colour was treated as a categorical variable, such that a hare on day d was either brown (pbrown), white (pwhite) or moulting (pmoult) and Σ(p1:3, j,d ) = 1. Camera site was coded as a random covariate si,j 2.5.2 | Temporal variation in moult phenology to reflect the hierarchical structure of the dataset and allow for repeated measurements. αi was the intercept and β1i was the effect of Next, to test the effect of annual variation in temperature and snow Julian day on the probability of being either brown, white or moult- season duration, we constructed an alternative set of univariate ing. Fall and spring moults were modelled separately. models. Each model included a single fixed effect of climate covariate β2i on the probability of being in a certain colour category to 2.5.1 | Spatial variation in moult phenology avoid multicollinearity issues (Supporting Information Table S1). The covariates included mean annual tmin, tmax, and duration of snow season in spring and fall at each camera site. The resulting β2i coef- First, to quantify average moult phenology of each population ficients were the slopes of reaction norms of the climate covariates (Canada, Colorado, New England), we combined colour observations on probabilities of being brown (β2brown) or white (β2white). from all years at that area and ran the model separately for each. We For all models, we obtained posterior distributions of all param- used the estimated probabilities to derive approximate dates when eters along with their 95% CRIs using Markov chain Monte Carlo �508 | ZIMOVA et al. (version 4.0.1), which we called using white-to-brown moults first, in late March, followed almost 2 weeks the R2jags package (Su & Yajima, 2012). Model convergence was later by populations in Canada and 4 weeks later by the southern- assessed using the Gelman–Rubin statistic, where values < 1.1 in- most population in Colorado. The Colorado population took the dicated convergence (Gelman & Rubin, 1992). We generated three shortest to complete the transitions (44 days versus 48 and 51 days chains of 300,000 iterations after a burn-in of 150,000 iterations in Canada and New England, respectively) and became brown only 2 and thinned by three. Parameters αi, β1i and β2i received a vague and 3 weeks later than the populations in Canada and New England, prior of N(0, 0.001), and the standard deviation of random effect si,j respectively (Figure 1b). (MCMC) implemented in jags received a uniform prior of U(0, 100). Variation in moult phenology between populations did not follow the north–south latitudinal gradient as predicted. Among popula- 2.5.3 | Camouflage mismatch tions, latitude had a significant effect on the spring moult phenology, but the effect was negative; hares at higher (i.e., northern) latitudes became brown earlier than hares in lower latitudes (Table 2). In the Camouflage mismatch was calculated based on the daily presence or fall, latitude had no effect on moult phenology (Table 2). Local cli- absence of snow and the modelled coat colour at each camera site. mate and elevation were strong predictors of moult phenology and Snow was present at a camera site when SWE > 0, and absent when SWE = 0, based on daily gridded 1 km × 1 km resolution data (SNODAS, for validation of dataset see Sirén et al., 2018). Next, we defined white hares as when mean pwhite ≥ 60% and brown hares as pbrown > 60% as these thresholds included mostly white or brown hares, respectively, when compared to observations. The camera days with brown and white hares were calculated using colour probabilities from the models that included the best annual climate predictor (effect sizes in Table 3) in order to account for inter-annual variation in phenology. To quantify the annual frequency of camouflage mismatch TA B L E 2 Effect of latitude, elevation and long-term climate covariates on snowshoe hare moult phenology. Mean effect sizes and 95% credible interval (CRI) estimates for slopes for univariate models including data from all years and populations combined. Betas indicate effects of covariates on the probability of the moult’s final colour category (β2brown in the spring, β2white in the fall). Snow is the duration of continuous snow season (days), tmax and tmin are the mean minimum and maximum temperature (°C) in springs and falls during 1980–2009. Asterisks indicate CRIs not overlapping 0. Values reflect standardized data within each population, we calculated the number of days at all Covariate Fall β2white Spring β2brown camera sites (camera days) when the colour of hares would either Latitude 0.566 0.689* (−0.136, 1.293) (0.376, 1.012) 2.165* −1.325* (1.450, 3.033) (−1.631, −1.039) 0.446 −0.809* (−0.214, 1.123) (−1.143, −0.492) −1.855* 0.776* (−2.479, −1.288) (0.440, 1.123) −2.370* 1.280* (−2.894, −1.909) (0.998, 1.579) match or mismatch against the background colour. ‘White mismatch’ occurred on days when hares were white and snow was Elevation absent at the site. ‘Brown mismatch’ occurred on days when hares were brown and snow was present (Mills et al., 2013). ‘Match’ occurred on days when hares were white (brown) and snow was present (absent). The proportion of white mismatch occurrence, for example, was calculated as the count of all camera days when hares would experience white mismatch, divided by the total number of camera days (i.e., number of camera sites at each area mul- Snow tmax tmin tiplied by the total number of days in a season). We calculated the mismatch occurrence for two 4-month periods when mismatch might occur at all three study areas; 1 February–31 May (spring) and 1 September–31 December (fall). All proportions were multiplied by 100 for interpretation in %. 3 | R E S U LT S 3.1 | Spatial variation in moult phenology Snowshoe hare moult phenology varied across the three study areas (Figure 1). In the fall, populations in Colorado and Canada initiated fall moults in early October, with some evidence that hares in Canada moulted faster (32 days total) than hares in Colorado (42 days; Figure 1b). The hare population in New England initiated fall moults almost 3 weeks later and took the longest to complete the moults (46 days). In the spring, hares in New England initiated the TA B L E 3 Effect of annual temperature and snow season duration on moult phenology in snowshoe hares. Betas are the slopes of reaction norms β2 (= mean effect size of annual climate covariate) and their 95% credible intervals (CRIs) on the probability of the moult’s final colour category. Snow is the duration of continuous snow season (days), tmax and tmin are the mean minimum and maximum temperature (°C) in spring and fall each year during 2010-2017. Asterisks indicate CRIs not overlapping zero. Values reflect standardized data Covariate Fall β2white Spring β2brown Snow annual 1.466* −1.627* (1.009, 1.929) (−1.969, −1.303) −2.070* 1.587* (−2.850, −1.432) (1.208, 2.003) −2.344* 1.273* (−2.943–1.845) (0.921, 1.655) tmax annual tmin annual �ZIMOVA et al. | 509 always in the predicted direction in both seasons; earlier fall and later than observed during 1980–2009 (Figure 3). Brown mismatch was rare spring moults were associated with sites that are generally colder, in the New England population, with the exception of fall 2016, with snowier and located at higher elevations. Elevation and long-term 6% brown mismatch (Supporting Information Figure S1). Snowshoe minimum temperature had the strongest effect on moult phenology hares in Colorado had lower proportions of white mismatch days than in both seasons (Table 2). in New England but reached 7% in the fall of 2016, which had a very short snow season (Figure 3). In the springs, brown mismatch was more 3.2 | Temporal variation in moult phenology common than white mismatch in Colorado (Table 4). In Canada, white We found evidence of population-level temperature- and snow-mediated plasticity in moult phenology. All annual temperature and snow covariates affected moult phenology in the predicted direction; moults occurred later in the spring and potentially earlier in the fall during colder and or snowier years (Table 3). In the spring, this annual variation in temperature and snow duration resulted in 2- to 3-week differences in mean population initiation and completion dates between some years in Canada and New England (Figure 2). In contrast, we found no significant differences between spring initiation and completion dates in the Colorado population. Furthermore, we did not detect any differences in the fall moult phenology start or end dates in any population (Supporting Information Figure S3). 3.3 | Camouflage mismatch The occurrence of camouflage mismatch varied between study areas and years (Figure 3, Supporting Information Figure S1, Table 4). White mismatch (white hare against snowless background) was relatively infrequent at all sites and occurred only during low snow years in New England and Colorado. Hares in New England experienced the highest frequency of white mismatch during both seasons (Table 4), with the highest proportions of camera days with white mismatch in fall 2015 (15%) and in spring 2016 (9%); seasons with 37–40 fewer snow days F I G U R E 2 Estimated annual spring moult initiation (i) and completion (c) dates in the studied hare populations in Canada, Colorado and New England. Points show mean date estimates and are coloured by the annual duration of spring snow season (in days). Horizontal lines show 95% credible intervals (CRIs; overlapping CRIs identify same dates) [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 3 Annual proportions of camera days with white mismatch occurrences plotted against anomalies in the number of snow days each season in Canada, Colorado and New England. Study area-specific anomaly in the number of snow days was calculated for each year as the difference between the mean number of snow days during each season and the mean number of snow days during 1980–2009 at all camera sites. Seasons when hares experienced very high white mismath are marked by the year inscriptions. Photo depicts white mismatched hare [Colour figure can be viewed at wileyonlinelibrary.com] �510 | ZIMOVA et al. TA B L E 4 Modelled mean percent of camera days with white and brown mismatch at each study area. Mean percent were calculated based on annual estimates (Canada: six falls and five springs, Colorado and New England: four springs and four falls). Standard deviations are given in parentheses Canada White Colorado Brown New England White Brown White Brown Spring 0.56 (0.55) 11.96 (4.31) 1.74 (1.34) 6.39 (0.68) 2.41 (4.36) 0.98 (0.73) Fall 0.55 (0.47) 9.54 (3.60) 2.69 (2.80) 1.11 (0.66) 4.72 (7.00) 2.42 (2.68) mismatch was very rare, exceeding 1% in only three out of the nine et al., 2014), mountain hare (Watson, 1963), stoats (Feder, 1990) seasons of monitoring (Figure 3), but brown mismatch was frequent and rock ptarmigan (Salomonsen, 1939), although all were examined (Table 4, Supporting Information Figure S1). over relatively small spatial scales or conclusions were based on opportunistic observations and low sample sizes. Likewise, for multiple 4 | D I S CU S S I O N colour moulting species, the global distribution of genetically determined winter brown and winter white coat colour morphs is driven by variation in snow cover duration (Mills et al., 2018). Using by-catch photographs from remote camera traps, we quan- We found no evidence that moult phenology variation was dis- tified the spatial and temporal variation in snowshoe hare moult tributed along a latitudinal (i.e., north–south) gradient. While the phenology and camouflage mismatch across nearly the full range effects of climate and elevation were strong in both seasons, the of environmental conditions experienced by the species. To our effects of latitude were not detected in the fall and in the opposite knowledge, this is the first study to evaluate these processes at such direction than expected in the spring (i.e., during the spring moult, resolution in any seasonally moulting species. First, local climate the probability of being brown increased with latitude, Table 2). We conditions were the main drivers of spatial and temporal variation in think that the positive effect of latitude in the spring was driven by moult phenology; hares in colder and snowier sites and in higher el- spatial variation in climate; the southernmost study area, Colorado, evations moult later in the spring and earlier in the fall. Second, hare was as cold as the northernmost area. Although our study lacked the populations exhibited temperature- and snow-mediated plasticity in spatial coverage to give a definite conclusion on the role of latitude moult phenology; we found strong evidence of later spring- and sug- in moult phenology variation, we showed that latitude alone can- gestive evidence of earlier fall-moults during colder and or snowier not predict moult phenology across a species range. Future studies years. Third, the occurrence of camouflage mismatch varied in space that include more replicates across the latitudinal gradient will help and time, but white mismatch was more common in areas charac- elucidate the interacting effects of latitude and local environmental terized by shallow, short-lasting snowpack. Finally, because areas conditions on seasonal moult phenology. with shallow snowpack are expected to face the steepest declines in snow cover duration (Lute, Abatzoglou, & Hegewisch, 2015; Ning & Bradley, 2015), we conclude that hares occupying those southern edge portions of their range are the most vulnerable to camouflage 4.2 | Phenotypic plasticity in response to temperature and snow mismatch under current and future climate change. We found that annual temperature and snow affect the moult phenol- 4.1 | Local climate drives variation in moult phenology ogy of snowshoe hares. However, this variation resulted in differences in population moult initiation and completion dates between years only in the spring (Figure 2). For example, hares in Canada became brown 21 days earlier in spring 2016, which was on average 2.9°C Our analysis based on 448 camera sites and spanning over 15 de- warmer (tmin) with a snow season 22 days shorter than spring 2014 grees of latitude showed that local climate influenced moult phe- (Figure 2 and Supporting Information Figure S2). Similar difference nology in snowshoe hares (Figure 1). Mean phenology in both was observed between the same 2 years in New England, where the seasons structured strongly along elevational, temperature, and 2.4°C higher tmin and 21-day shorter snow season in spring 2016 snow cover gradients. Similar findings were described for other than in 2014 corresponded with 11- and 14-day advances in moult phenological traits across taxa (e.g., migration, reproduction and initiation and completion dates, respectively (Figure 2 and Supporting hibernation) where phenology correlated with local variation in cli- Information Figure S2). In contrast, we did not detect any differences mate (e.g., Duursma, Gallagher, & Griffith, 2018; Fielding, Whittaker, between the spring moult phenologies in the Colorado population Butterfield, & Coulson, 1999; Sheriff et al., 2011; While & Uller, during all years of monitoring (Figure 2). However, annual temperature 2014). As for seasonal colour moult phenology, local climatic factors and snow season length were less variable in Colorado than in the have been previously described to determine phenology in colour other two study areas (Supporting Information Figure S2). For exam- moulting species including snowshoe hares (Grange, 1932; Zimova ple, the 11-day difference in snow duration and 1.3°C difference in �| ZIMOVA et al. 511 tmin between the two most extreme springs (2011 and 2015) might for individual snowshoe hares in the Rocky Mountains in Montana not have been sufficient to observe significant differences between (i.e., 2–3 weeks in the spring; Mills et al., 2013; Zimova et al., 2014). phenology dates (Figure 2). To determine whether hares in Colorado However, in this study, the number of independent snowshoe hare have lower phenotypic plasticity or whether this finding was caused colour observations was lower during the falls than in the springs at by lower inter-annual variation will require analyses of additional years most of the sites (Table 1), which may have resulted in reduced statis- across a wider range of climatic conditions. tical power to detect differences during this season. We found some evidence for population-level phenotypic plas- Finally, we note that the methodology prevented phenology ticity in the fall moult phenology (Table 3), but we did not detect monitoring of the same individuals – and their fates – over multiple significant differences between initiation or completion dates in any years. Therefore, we were unable to investigate to what extent nat- study area (Supporting Information Figure S3). Lower phenological ural selection against camouflage mismatch might have contributed plasticity in the fall than in the spring was previously described in to the observed differences in temporal variation moult phenology. least weasels and connected to lower variation in the fall than in the Specifically, during the few very short snow duration years, natural spring temperatures (Atmeh et al., 2018). In this study, the lower in- selection might have removed highly mismatched individuals from the ter-annual variation in climatic conditions might have similarly con- population and contributed to the observed shift in the phenological tributed to the seasonal variation in plasticity in some cases (e.g., low distribution. Quantifying the relative importance of natural selection variation in tmax in fall in Canada; Supporting Information Figure and individual plasticity on population-level responses in phenology S2). However, in most cases, the inter-annual variation in tempera- will require intensive field studies with tagged individuals over mul- ture and snow at our study areas was comparable between spring tiple years. and fall seasons (Supporting Information Figure S2). For example, snow season duration differed by up to 44 days and tmax by 2.3°C between the most extreme falls in New England (2015 and 2016; 4.3 | Spatial variation in camouflage mismatch Supporting Information Figure S2), yet we observed no differences between moult phenology dates (Supporting Information Figure S3). We observed high variation in the number of snow days each year but, Furthermore, low plasticity in the fall and a similar level of phenotypic overall, snow cover duration has decreased over the last 30 years in plasticity in the spring as observed here were previously described all three areas (Supporting Information Figure S4; Fassnacht, Venable, F I G U R E 4 Mean daily snow water equivalent (SWE; mm) at the remote camera sites for the years of moult phenology monitoring. Coloured circles along the x axes indicate the population mean moult completion dates for each year, with spring moults on the left and fall moults on the right. Mean number of snow days for years of monitoring is given in facet titles for spring and fall [Colour figure can be viewed at wileyonlinelibrary.com] �512 | ZIMOVA et al. McGrath, & Patterson, 2018; Harpold et al., 2012; Kunkel et al., 2016; when hares were white than vice versa (Figure 4). In contrast, during Mote, Li, Lettenmaier, Xiao, & Engel, 2018; Ning & Bradley, 2015). The springs in Canada and Colorado, hares experienced long periods declines in snow cover manifested differently in each population, how- when they had already moulted to summer brown pelage, but snow ever. In the Canadian population, white mismatch was very low each was still present on the ground (Figure 4). Furthermore, the onset year, while the number of snow days ranged from 1 to 28 fewer days of snow during fall often occurred prior to hares completing their than observed on average during 1980–2009 (Figure 3). This is likely brown-to-white moults in Canada. Overall, the relatively high occur- due to the very deep, long-lasting snowpack in the study area (Table 1, rence of brown camouflage mismatch, despite the recorded snow Figure 4). For example, during the years of monitoring, spring snow declines, was unexpected. Eventually, however, as climate change season ended on average on 9 June (Figure 4), more than 2 weeks will continue to lead to shorter snow duration across most of the after hares finished their white-to-brown moults (i.e., moults initiated snowshoe hare range (Danco, DeAngelis, Raney, & Broccoli, 2016; on 5 April and completed on 23 May; Figure 1b). Yet, annual snow Easterling et al., 2017; Fyfe et al., 2017), the occurrence of brown season duration in the Canadian Rockies is predicted to decrease by mismatch will decrease. about 1 month by the end of the century (Pomeroy, Fang, & Rasouli, The relative fitness costs of white versus brown mismatch are 2015). Those predictions would advance average snowmelt to about unknown, but we suspect that white mismatch has a higher survival 2 weeks prior to when hares become fully brown. Although increase in cost than brown mismatch based on our experience in the field, white camouflage mismatch will likely occur in very short snow years, while locating radio-collared hares. First, brown animals and objects based on our results a mean decrease of snow season by 1 month will (e.g., branches, tree trunks, brown animals) are relatively common unlikely lead to dramatic increases in the frequency of white hares on year-round, but white animals and objects are rare outside of win- snowless background (see Figure 3). ter. Perhaps due in part to this frequency difference in the two mis- In contrast, hares in Colorado and New England are much more match types, a white hare against a snowless background appears vulnerable to white camouflage mismatch. For both populations, the far more conspicuous than a brown hare resting on snow. Previous proportion of white mismatch increased rapidly with fewer snow quantifications of survival costs were carried out for ‘absolute mis- days, especially once the snow days anomaly exceeded 21 days match’, that is both white and brown mismatch combined (Wilson (Figure 3). Beyond this 3-week threshold, hares began to experience et al., 2018; Zimova et al., 2016). Nonetheless, as documented here elevated white mismatch (i.e., 15% in fall and 9% in spring in New and elsewhere (Wilson et al., 2018; Zimova et al., 2016), white mis- England; 7% in fall in Colorado). This suggests that phenotypic plas- match is already high in some populations and will increase under cli- ticity in moult phenology is insufficient to buffer against the snow mate change (Mills et al., 2013). Therefore, fitness costs of white and declines in those marginal areas where snow season is already short brown mismatch should be quantified to inform conservation efforts, (Table 1, Figure 4). To contrast the New England study area with notably in situ management actions that foster evolutionary rescue, the previous example from Canada, the continuous snow season or genetic rescue by assisted gene flow of individuals with pre- ended on average on 18 April (Figure 4) but hares underwent spring adapted moult phenologies or winter coat colour (Mills et al., 2018). moults from 23 March to 13 May (Figure 1b). Therefore, with mean Understanding the spatial and temporal variation in phenologi- snowmelt occurring before hares are halfway through the moult, cal traits is critical for understanding the impact of climate change early snowmelt years (e.g., 2016 seasonal snow melted on 1 April; and species’ adaptive potential to environmental stressors. Here, Supporting Information Figure S2) result in steep increases in white we showed that snowshoe hare moult phenology is determined by mismatch. local climate, but populations vary in their susceptibility to camou- Based on our results, in the absence of evolutionary shifts in flage mismatch. Snowshoe hares responded to annual variation in phenology, snowshoe hares in New England and Colorado will temperature and snow via some adjustments in moult phenology, suffer drastic increases in white mismatch over the next century. but the buffering effects of phenotypic plasticity were diminished In the Colorado study area, the frequency of low snow years is in populations distributed along the southern edge of their range. In expected to increase and the annual snow season to decrease by those areas, characterized by mild climate and shallow, short-lived 20–50 days (for comparison see mismatch in fall 2016, Figure 3) snowpack, climate change mediated snow declines led to higher by mid-century (Lute et al., 2015). Similarly, in many parts of New phenological mismatch. Therefore, hares occupying southern, mar- England, snow cover duration is predicted to shorten from the ginal areas will in the absence of rapid evolution experience steep current 5 months to 3 months by 2100 (Ning & Bradley, 2015). increases in camouflage mismatch, as those areas are expected to Under such projections, hare in both study areas will by the end experience the largest declines in snow cover duration (Fyfe et al., of the century routinely experience the same amount of mismatch 2017; Ning & Bradley, 2015), consistent with theoretical expecta- as they experienced during the very low snow years of fall 2016 in tions of range contraction (Sirén & Morelli, 2019). More generally, Colorado and 2015 in New England. our results underscore that populations vary in their susceptibility The pattern observed with white mismatch was somewhat mir- to environmental stressors and management efforts should consider rored by that of brown mismatch. In both seasons in New England this intraspecific variation to identify populations most vulnerable and during the fall in Colorado, white mismatch was more frequent under global environmental change (Hampe & Petit, 2005; Nadeau, than brown mismatch as snow cover was more likely to be absent Urban, & Bridle, 2017). �| ZIMOVA et al. AC K N OW L E D G M E N T S We are grateful to the many field technicians and volunteers who managed the remote camera sites and photographs. Next, we thank Sean T. Giery, Zac A. Cheviron and two anonymous referees for their comments on earlier drafts. This work was supported by the Department of the Interior Southeast Climate Adaptation Science Center Global Change Fellowship through Cooperative Agreement No. G10AC00624 to MZ; Department of the Interior Northeast Climate Adaptation Science Center, T-2-3R grant for Nongame Species Monitoring and Management through the New Hampshire Fish and Game Department, and a E-1-25 grant for Investigations and Population Recovery through the Vermont Fish and Wildlife Department to APKS; Colorado Species Conservation Trust Fund; the National Science Foundation Division of Environmental Biology Grant 0841884 and 1907022 to LSM; the National Science Foundation EPSCoR Award No. 1736249 (OIA-1736249); North Carolina State University; and the University of Montana. Parks Canada funded data collection in Banff, Kootenay and Yoho National Parks. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. DATA ACC E S S I B I L I T Y The data used in this study are available at datadryad.org (https​:// doi.org/10.5061/dryad.k98sf​7m35). ORCID https://orcid.org/0000-0002-8264-9879 Marketa Zimova Alexej P. K. Sirén https://orcid.org/0000-0003-3067-6418 Joshua J. Nowak https://orcid.org/0000-0002-8553-7450 Jacob S. Ivan https://orcid.org/0000-0002-5314-7757 https://orcid.org/0000-0001-5865-5294 Toni Lyn Morelli Jesse Whittington L. Scott Mills https://orcid.org/0000-0002-4129-7491 https://orcid.org/0000-0001-8771-509X REFERENCES Atmeh, K., Andruszkiewicz, A., & Zub, K. (2018). Climate change is affecting mortality of weasels due to camouflage mismatch. Scientific Reports, 8, 7648. https​://doi.org/10.1038/s41598-018-26057-5 Barrett, A. P. (2003). National operational hydrologic remote sensing center snow data assimilation system (SNODAS) products at NSIDC. Boulder, CO: National Snow and Ice Data Center. Both, C., Bouwhuis, S., Lessells, C. M., & Visser, M. E. (2006). Climate change and population declines in a long-distance migratory bird. Nature, 441, 81–83. https​://doi.org/10.1038/natur​e 04539 Bradshaw, W. E., & Holzapfel, C. M. (2007). Evolution of animal photoperiodism. Annual Review of Ecology Evolution and Systematics, 38, 1–25. https​://doi.org/10.1146/annur​ev.ecols​ys.37.091305.110115 Cheng, E., Hodges, K. E., Melo-Ferreira, J., Alves, P. C., & Mills, L. S. (2014). Conservation implications of the evolutionary history and genetic diversity hotspots of the snowshoe hare. Molecular Ecology, 23, 2929–2942. https​://doi.org/10.1111/mec.12790​ Choi, G., Robinson, D. A., & Kang, S. (2010). Changing Northern Hemisphere snow seasons. Journal of Climate, 23, 5305–5310. https​://doi.org/10.1175/2010J​CLI36​4 4.1 Chuine, I. (2010). Why does phenology drive species distribution? Philosophical Transactions of the Royal Society B: Biological Sciences, 365, 3149–3160. https​://doi.org/10.1098/rstb.2010.0142 513 Cohen, J. M., Lajeunesse, M. J., & Rohr, J. R. (2018). A global synthesis of animal phenological responses to climate change. Nature Climate Change, 8, 224–228. https​://doi.org/10.1038/s41558-018-0067-3 Danco, J. F., DeAngelis, A. M., Raney, B. K., & Broccoli, A. J. (2016). Effects of a warming climate on daily snowfall events in the Northern Hemisphere. Journal of Climate, 29, 6295–6318. https​:// doi.org/10.1175/JCLI-D-15-0687.1 Diefenbach, D. R., Rathbun, S. L., Vreeland, J. K., Grove, D., & Kanapaux, W. J. (2016). Evidence for range contraction of snowshoe hare in Pennsylvania. Northeastern Naturalist, 23, 229–248. https​://doi. org/10.1656/045.023.0205 Dietz, A. J., Kuenzer, C., & Dech, S. (2015). Global SnowPack: a new set of snow cover parameters for studying status and dynamics of the planetary snow cover extent. Remote Sensing Letters, 6, 844–853. Duursma, D. E., Gallagher, R. V., & Griffith, S. C. (2018). Variation in the timing of avian egg-laying in relation to climate. Ecography, 42, 535–548. Easterling, D. R., Kunkel, K. E., Arnold, J. R., Knutson, T., LeGrande, A. N., Leung, L. R., … Wehner, M. F. (2017). Precipitation change in the United States. In D. J. Wuebbles, D. W. Fahey, K. A. Hibbard, D. J. Dokken, B. C. Stewart, & T. K. Maycock (Eds.), Climate Science Special Report: Fourth National Climate Assessment (Vol. 1, pp. 207–230. U.S. Global Change Research Program. Fassnacht, S. R., Venable, N. B. H., McGrath, D., & Patterson, G. G. (2018). Sub-seasonal snowpack trends in the Rocky Mountain National Park Area, Colorado, USA. Water, 10, 1–19. https​://doi.org/10.3390/ w1005​0562 Feder, S. (1990). Environmental determinants of seasonal coat color change in weasels (Mustela erminea) from two populations (MS Thesis). University of Alaska Fairbanks. Fielding, C. A., Whittaker, J. B., Butterfield, J. E. L., & Coulson, J. C. (1999). Predicting responses to climate change: The effect of altitude and latitude on the phenology of the Spittlebug Neophilaenus lineatus. Functional Ecology, 13, 65–73. https​://doi. org/10.1046/j.1365-2435.1999.00009.x Fyfe, J. C., Derksen, C., Mudryk, L., Flato, G. M., Santer, B. D., Swart, N. C., … Jiao, Y. (2017). Large near-term projected snowpack loss over the western United States. Nature Communications, 8, 14996. https​:// doi.org/10.1038/ncomm​s14996 Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–511. https​://doi. org/10.1214/ss/11770​11136​ Gienapp, P., & Brommer, J. E. (2014). Evolutionary dynamics in response to climate change. In A. Charmantier, D. Garant, & L. E. B. Kruuk (Eds.), Quantitative genetics in the wild (pp. 254–273). Oxford, UK: Oxford University Press. Gigliotti, L. C., Diefenbach, D. R., & Sheriff, M. J. (2017). Geographic variation in winter adaptations of snowshoe hares (Lepus americanus). Canadian Journal of Zoology, 95, 539–545. Grange, W. B. (1932). The pelages and color changes of the snowshoe hare, Lepus americanus phaeonotus, Allen. Journal of Mammalogy, 13, 99–116. https​://doi.org/10.2307/1374046 Hampe, A., & Petit, R. J. (2005). Conserving biodiversity under climate change: The rear edge matters. Ecology Letters, 8, 461–467. https​:// doi.org/10.1111/j.1461-0248.2005.00739.x Harpold, A., Brooks, P., Rajagopal, S., Heidbuchel, I., Jardine, A., & Stielstra, C. (2012). Changes in snowpack accumulation and ablation in the intermountain west. Water Resources Research, 48, W11501. https​://doi.org/10.1029/2012W​R011949 Hodges, K. E. (1999). Proximate factors affecting snowshoe hare movements during a cyclic population low phase. Ecoscience, 6, 487–496. https​://doi.org/10.1080/11956​860.1999.11682558 Hofman, M. A. (2004). The brain’s calendar: Neural mechanisms of seasonal timing. Biological Reviews, 79, 61–77. https​://doi.org/10.1017/ S1464​79310​3 006250 �514 | Hut, R. A., Paolucci, S., Dor, R., Kyriacou, C. P., & Daan, S. (2013). Latitudinal clines: An evolutionary view on biological rhythms. Proceedings of the Royal Society B: Biological Sciences, 280, 20130433. Imperio, S., Bionda, R., Viterbi, R., & Provenzale, A. (2013). Climate change and human disturbance can lead to local extinction of Alpine rock ptarmigan: New insight from the western Italian Alps. PLoS ONE, 8, e81598. https​://doi.org/10.1371/journ​al.pone.0081598 Ivan, J. S. (2011). Density, demography, and seasonal movements of snowshoe hares in central Colorado (PhD Dissertation). Colorado State University. Krebs, C. J., Boonstra, R., Boutin, S., & Sinclair, A. R. E. (2001). What drives the 10-year cycle of snowshoe hares? BioScience, 51, 25–35. https​:// doi.org/10.1641/0006-3568(2001)051[0025:WDTYC​O]2.0.CO;2 Kunkel, K. E., Robinson, D. A., Champion, S., Yin, X., Estilow, T., & Frankson, R. M. (2016). Trends and extremes in Northern Hemisphere snow characteristics. Current Climate Change Reports, 2, 65–73. https​://doi.org/10.1007/s40641-016-0036-8 Lane, J. E., Kruuk, L. E. B., Charmantier, A., Murie, J. O., & Dobson, F. S. (2012). Delayed phenology and reduced fitness associated with climate change in a wild hibernator. Nature, 489, 554–557. https​://doi. org/10.1038/natur​e11335 Livneh, B., Bohn, T. J., Pierce, D. W., Munoz-Arriola, F., Nijssen, B., Vose, R., … Brekke, L. (2015). A spatially comprehensive, hydrometeorological data set for Mexico, the U.S., and Southern Canada 1950–2013. Scientific Data, 2, 150042. https​://doi.org/10.1038/ sdata.2015.42 Lute, A. C., Abatzoglou, J. T., & Hegewisch, K. C. (2015). Projected changes in snowfall extremes and interannual variability of snowfall in the western United States. Water Resources Research, 51, 960–972. https​://doi.org/10.1002/2014W​R016267 Mills, L. S., Bragina, E. V., Kumar, A. V., Zimova, M., Lafferty, D. J. R., Feltner, J., … Fay, K. (2018). Winter color polymorphisms identify global hot spots for evolutionary rescue from climate change. Science, 359, 1033–1036. https​://doi.org/10.1126/scien​ce.aan8097 Mills, L. S., Griffin, P. C., Hodges, K. E., McKelvey, K., Ruggiero, L., & Ulizio, T. (2005). Pellet count indices compared to mark-recapture estimates for evaluating snowshoe hare density. Journal of Wildlife Management, 69, 1053–1062. https​://doi.org/10.2193/0022-5 41X(2005)069[1053:PCICT​M]2.0.CO;2 Mills, L. S., Zimova, M., Oyler, J., Running, S., Abatzoglou, J. T., & Lukacs, P. M. (2013). Camouflage mismatch in seasonal coat color due to decreased snow duration. Proceedings of the National Academy of Sciences USA, 110, 7360–7365. https​://doi.org/10.1073/pnas.12227​24110​ Mote, P. W., Li, S., Lettenmaier, D. P., Xiao, M., & Engel, R. (2018). Dramatic declines in snowpack in the western US. Climate and Atmospheric Science, 1, 2. https​://doi.org/10.1038/s41612-018-0012-1 Nadeau, C. P., Urban, M. C., & Bridle, J. R. (2017). Climates past, present, and yet-to-come shape climate change vulnerabilities. Trends in Ecology and Evolution, 32, 786–800. https​://doi.org/10.1016/ j.tree.2017.07.012 Nagorsen, D. W. (1983). Winter pelage colour in snowshoe hares (Lepus americanus) from the Pacific Northwest. Canadian Journal of Zoology, 61, 2313–2318. Ning, L., & Bradley, R. S. (2015). Snow occurrence changes over the central and eastern United States under future warming scenarios. Scientific Reports, 5, 17073. https​://doi.org/10.1038/srep1​7073 Parmesan, C., & Yohe, G. (2003). A globally coherent fingerprint of climate change impacts across natural systems. Nature, 421, 37–42. https​://doi.org/10.1038/natur​e 01286 Pedersen, S., Odden, M., & Pedersen, H. C. (2017). Climate change induced molting mismatch? Mountain hare abundance reduced by duration of snow cover and predator abundance. Ecosphere, 8, e01722. https​://doi.org/10.1002/ecs2.1722 ZIMOVA et al. Phillimore, A. B., Hadfield, J. D., Jones, O. R., & Smithers, R. J. (2010). Differences in spawning date between populations of common frog reveal local adaptation. Proceedings of the National Academy of Sciences USA, 107, 8292–8297. Pomeroy, J. W., Fang, X., & Rasouli, K. (2015). Sensitivity of snow processes to warming in the Canadian Rockies. In Proceedings of the 72nd Eastern Snow Conference, 9–11 June 2015, (pp. 22–33). Sherbrooke, Quebec, Canada. Post, E., & Forchhammer, M. C. (2008). Climate change reduces reproductive success of an Arctic herbivore through trophic mismatch. Philosophical Transactions of the Royal Society B: Biological Sciences, 363, 2369–2375. https​://doi.org/10.1098/rstb.2007.2207 R Core Team (2016). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from http://www.R-proje​c t.org/. Salomonsen, F. (1939). Moults and sequence of plumages in the rock ptarmigan (Lagopus mutus (Montin)). Videnskabelige Meddelelser, 103, 1–491. Senner, N. R., Stager, M., & Sandercock, B. K. (2017). Ecological mismatches are moderated by local conditions for two populations of a long-distance migratory bird. Oikos, 126, 61–72. https​://doi. org/10.1111/oik.03325​ Sheriff, M. J., Kenagy, G. J., Richter, M., Lee, T., Toien, O., Kohl, F., … Barnes, B. M. (2011). Phenological variation in annual timing of hibernation and breeding in nearby populations of Arctic ground squirrels. Proceedings of the Royal Society B: Biological Sciences, 278, 2369–2375. https​://doi.org/10.1098/rspb.2010.2482 Sirén, A. P. K., & Morelli, T. L. (2019). Interactive range-limit theory (iRLT): A theoretical extension for predicting ranges shifts. Journal of Animal Ecology, 1–15. https​://doi.org/10.1111/1365-2656.13150​ Sirén, A. P. K., Somos-Valenzuela, M., Callahan, C., Kilborn, J. R., Duclos, T., Tragert, C., & Morelli, T. L. (2018). Looking beyond wildlife: Using remote cameras to evaluate accuracy of gridded snow data. Remote Sensing in Ecology and Conservation, 4, 375–386. https​://doi. org/10.1002/rse2.85 Su, Y. S., & Yajima, M. (2012). R2jags: A package for running jags from R. Retrieved from http://CRAN.R-proje​c t.org/packa​ge=R2jags. Sultaire, S. M., Pauli, J. N., Martin, K. J., Meyer, M. W., Notaro, M., & Zuckerberg, B. (2016). Climate change surpasses land-use change in the contracting range boundary of a winter-adapted mammal. Proceedings of the Royal Society of London B: Biological Sciences, 283, 20153104. Thackeray, S. J., Sparks, T. H., Frederiksen, M., Burthe, S., Bacon, P. J., Bell, J. R., … Wanless, S. (2010). Trophic level asynchrony in rates of phenological change for marine, freshwater and terrestrial environments. Global Change Biology, 16, 3304–3313. https​://doi. org/10.1111/j.1365-2486.2010.02165.x Thornton, P. E., Thornton, M. M., Mayer, B. W., Wei, Y., Devarakonda, R., Vose, R. S., & Cook, R. B. (2018). Daymet: Daily surface weather data on a 1-km Grid for North America, version 3. ORNL DAAC, Oak Ridge, TN. Retrieved from https​://doi.org/10.3334/ORNLD​ AAC/1328. Vaughan, D. G., Comiso, J., Allison, I., Carrasco, J., Kaser, G., Kwok, R., & Holland, D. (2013). Observations: Cryosphere. Climate change 2013: The physical science basis. Contribution of working group I to the fifth assessment report of the Intergovernmental Panel on Climate Change. UK: Cambridge University Press. Visser, M. E., Caro, S. P., van Oers, K., Schaper, S. V., & Helm, B. (2010). Phenology, seasonal timing and circannual rhythms: Towards a unified framework. Philosophical Transactions of the Royal Society B: Biological Sciences, 365, 3113–3127. Visser, M. E., & Gienapp, P. (2019). Evolutionary and demographic consequences of phenological mismatches. Nature Ecology and Evolution, 3, 879–885. https​://doi.org/10.1038/s41559-019-0880-8 �| ZIMOVA et al. Watson, A. (1963). The effect of climate on the colour changes of mountain hares in Scotland. Proceedings of the Zoological Society of London, 141, 823–835. While, G. M., & Uller, T. (2014). Quo vadis amphibia? Global warming and breeding phenology in frogs, toads and salamanders. Ecography, 37, 921–929. https​://doi.org/10.1111/ecog.00521​ Williams, C. M., Henry, H. A. L., & Sinclair, B. J. (2015). Cold truths: How winter drives responses of terrestrial organisms to climate change. Biological Reviews, 90, 214–235. https​://doi.org/10.1111/brv.12105​ Wilson, E. C., Shipley, A. A., Zuckerberg, B., Peery, M. Z., & Pauli, J. N. (2018). An experimental translocation identifies habitat features that buffer camouflage mismatch in snowshoe hares. Conservation Letters, 12, e12614. https​://doi.org/10.1111/conl.12614​ Zimova, M., Barnard, L., Davis, B. M., Kumar, A. V., Lafferty, D. J. R., & Mills, L. S. In Review. Using remote cameras to measure seasonal molts. Ecosphere. Zimova, M., Hackländer, K., Good, J. M., Melo-Ferreira, J., Alves, P. C., & Mills, L. S. (2018). Function and underlying mechanisms of seasonal colour moulting in mammals and birds: What keeps them changing in a warming world? Biological Reviews, 93, 1478–1498. https​://doi. org/10.1111/brv.12405​ Zimova, M., Mills, L. S., Lukacs, P. M., & Mitchell, M. S. (2014). Snowshoe hares display limited phenotypic plasticity to mismatch in seasonal camouflage. Proceedings of the Royal Society B: Biological Sciences, 281, 20140029. https​://doi.org/10.1098/rspb.2014.0029 Zimova, M., Mills, L. S., & Nowak, J. J. (2016). High fitness costs of climate change-induced camouflage mismatch. Ecology Letters, 19, 299–307. https​://doi.org/10.1111/ele.12568​ 515 B I O S K E TC H Marketa Zimova is a postdoctoral researcher interested in the effects of global environmental change on wild populations. The majority of her work aims to understand species adaptive potential to climate change and to yield useful information that can guide conservation and management decisions necessary for maintaining biodiversity. The research team includes scientists from the United States and Canada who study climate change and the ecology and evolutionary biology of wild organisms. S U P P O R T I N G I N FO R M AT I O N Additional supporting information may be found online in the Supporting Information section. How to cite this article: Zimova M, Sirén APK, Nowak JJ, et al. Local climate determines vulnerability to camouflage mismatch in snowshoe hares. Global Ecol Biogeogr. 2020;29:503–515. https​://doi.org/10.1111/geb.13049​ � https://cpw.cvlcollections.org/files/original/a9aea2b9b51eb41875785c0393a2cd97.pdf f949c23e5bd748968865f9449de16480 PDF Text Text Supporting Information for: LOCAL CLIMATE DETERMINES VULNERABILITY TO CAMOUFLAGE MISMATCH IN SNOWSHOE HARES Marketa Zimova, Alexej Sirén, Joshua Nowak, Alexander M. Bryan, Jacob S. Ivan, Toni Lyn Morelli, J. Skyler Suhrer, Jesse Whittington, L. Scott Mills Figures S1-S4 Tables S1-S3 1 �Figure S1. Annual proportions of camera days with white (red) and brown (gray) mismatch occurrences (in %) plotted against anomalies in the number of snowdays each season in Canada, Colorado and New England. Study area-specific anomaly in the number of snow days was calculated for each year as the difference between the mean number of snow days during each season and the mean number of snow days during 1980-2009 at all camera sites. 2 �a b c Figure S2. Annual mean and standard deviation in snow and temperature variables at the three study areas during the years of monitoring and for the 30-year mean (black points; 1980-2009). (a) duration of snow season (days; falls = bottom; spring = top). (b) seasonal minimum (bottom) and maximum (top) temperatures in spring. (c) seasonal minimum (bottom) and maximum (top) temperature in fall. Annotations at the bottom right corner of each facet give standard deviation in the depicted climate variable across the years of monitoring. 3 �Figure S3. Estimated annual fall molt initiation (left) and completion (right) dates in the hare populations in Canada, Colorado and New England. Points show mean date estimates and are colored by the annual duration of fall snow season (in days). Horizontal lines show 95% credible intervals (overlapping CRIs identify same dates). 4 �Figure S4. Mean annual number of days with snow at camera trap sites in the fall (bottom series) and spring (top series) 1980-2017. The colored points highlight years of molt phenology monitoring. Vertical lines show standard deviations. Dashed lines show linear regression slopes for each study area. Numbers in the bottom show estimated reductions in number of snow days between 1980-2017. 5 �Table S1. Pearson’s correlation coefficients between annual minimum (tmin) or maximum (tmax) temperature and seasonal snow duration during the course of the study (2010-2017) and the 30-year period (1980-2009). The correlation coefficients are calculated across all camera sites. Fall Spring Time Period tmin x tmax tmin x snow tmax x snow tmin x tmax tmin x snow tmax x snow 2010-2017 1980-2009 0.77 0.66 -0.68 -0.85 -0.73 -0.83 0.78 0.61 -0.62 -0.74 -0.66 -0.64 6 �Table S2. Effects of geospatial and long-term climate covariates on snowshoe hare molt phenology. Mean effect sizes and 95% credible intervals (CRI) estimates for slopes for univariate models including data from all years and populations combined. Betas indicate effect of latitude, elevation, duration of snow season, mean seasonal minimum (tmin) and maximum (tmax) temperature during 1980-2009 on the probability of brown (b2brown) and white (b2white) coat color. Asterisks indicate CRIs not overlapping 0. Values reflect standardized data. Covariate Fall b2brown Latitude Elevation Snow tmax tmin Fall b2white Spring b2brown Spring b2white -0.624 0.566 0.689* -1.127* (-1.474, 0.165) (-0.136, 1.293) (0.376, 1.012) (-1.500, -0.764) -1.952* 2.165* -1.325* 1.812* (-2.739, -1.224) (1.450, 3.033) (-1.631, -1.039) (1.513, 2.140) -0.752* -0.809* 0.382* 0.446 (-1.494, -0.077) (-0.214, 1.123) ( -1.143, -0.492) (0.015, 0.746) 1.693* -1.855* 0.776* (1.046, 2.409) (-2.479, -1.288) (0.440, 1.123) (-0.938, -0.165) 2.182* -2.370* -1.665* (1.542, 2.976) (-2.894, -1.909) (0.998, 1.579) 1.280* -0.546* (-2.029, -1.325) 7 �Table S3. Effects of annually varying covariates on snowshoe hare molt phenology. Betas are the slopes of reaction norms b2 (=mean effect size of annually varying climate covariate) and their 95% credible intervals (CRI) on probability of being in brown (b2brown) and white coat color (b2white). Results are based on univariate models using standardized data. Asterisks indicate CRIs not overlapping zero. Covariate Snow annual tmax annual tmin annual Fall b2brown Fall b2white Spring b2brown Spring b2white -1.036* 1.466* -1.627* 1.965* (-1.618, -0.475) (1.009, 1.929) (-1.969, -1.303) (1.637, 2.320) 1.718* -2.070* 1.587* -0.831* (1.063, 2.489) (-2.850, -1.432) (1.208, 2.003) (-1.239, -0.434) 2.122* -2.344* 1.273* -1.090* (1.427, 2.944) (-2.943 -1.845) (0.921, 1.655) (-1.480, -0.718) 8 � Dublin Core The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/. Title A name given to the resource Journal Articles Description An account of the resource CPW peer-reviewed journal publications Text A resource consisting primarily of words for reading. Examples include books, letters, dissertations, poems, newspapers, articles, archives of mailing lists. Note that facsimiles or images of texts are still of the genre Text. Dublin Core The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/. Title A name given to the resource Local climate determines vulnerability to camouflage mismatch in snowshoe hares Description An account of the resource <strong>Aim</strong> <br />Phenological mismatches, when life-events become mistimed with optimal environmental conditions, have become increasingly common under climate change. Population-level susceptibility to mismatches depends on how phenology and phenotypic plasticity vary across a species’ distributional range. Here, we quantify the environmental drivers of colour moult phenology, phenotypic plasticity, and the extent of phenological mismatch in seasonal camouflage to assess vulnerability to mismatch in a common North American mammal.<br /> <p class="article-section__sub-title section1"><strong>Location</strong><br />North America.</p> <p class="article-section__sub-title section1"><strong>Time period</strong><br />2010–2017.</p> <p class="article-section__sub-title section1"><strong>Major taxa studied</strong><br />Snowshoe hare (<i>Lepus americanus</i>).</p> <p class="article-section__sub-title section1"><strong>Methods</strong><br />We used &gt; 5,500 by-catch photographs of snowshoe hares from 448 remote camera trap sites at three independent study areas. To quantify moult phenology and phenotypic plasticity, we used multinomial logistic regression models that incorporated geospatial and high-resolution climate data. We estimated occurrence of camouflage mismatch between hares’ coat colour and the presence and absence of snow over 7 years of monitoring.</p> <p class="article-section__sub-title section1"><strong>Results</strong><br />Spatial and temporal variation in moult phenology depended on local climate conditions more so than on latitude. First, hares in colder, snowier areas moulted earlier in the fall and later in the spring. Next, hares exhibited phenotypic plasticity in moult phenology in response to annual variation in temperature and snow duration, especially in the spring. Finally, the occurrence of camouflage mismatch varied in space and time; white hares on dark, snowless background occurred primarily during low-snow years in regions characterized by shallow, short-lasting snowpack.</p> <p class="article-section__sub-title section1"><strong>Main conclusions</strong><br />Long-term climate and annual variation in snow and temperature determine coat colour moult phenology in snowshoe hares. In most areas, climate change leads to shorter snow seasons, but the occurrence of camouflage mismatch varies across the species’ range. Our results underscore the population-specific susceptibility to climate change-induced stressors and the necessity to understand this variation to prioritize the populations most vulnerable under global environmental change.</p> Bibliographic Citation A bibliographic reference for the resource. Recommended practice is to include sufficient bibliographic detail to identify the resource as unambiguously as possible. Zimova, M., A. P. Sirén, J. J. Nowak, A. M. Bryan, J. S. Ivan, T. L. Morelli, S. L. Suhrer, J. Whittington, and L. S. Mills. 2019. Local climate determines vulnerability to camouflage mismatch in snowshoe hares. Global Ecology and Biogeography 29:503–515.  <a href="https://doi.org/10.1111/geb.13049" target="_blank" rel="noreferrer noopener">https://doi.org/10.1111/geb.13049</a> Creator An entity primarily responsible for making the resource Zimova, Marketa Sirén, Alexej P. K. Nowak, Joshua J. Bryan, Alexander M. Ivan, Jacob S. Morelli, Toni Lyn Suhrer, Skyler L. Whittington, Jesse Mills, L. Scott Subject The topic of the resource Adaptation Camouflage mismatch Climate change Latitudinal gradient Phenological mismatch Phenotypic plasticity Range edge Snow Snowshoe hares Extent The size or duration of the resource. 13 pages Date Created Date of creation of the resource. 2019-12-26 Rights Information about rights held in and over the resource <a href="http://rightsstatements.org/vocab/InC-NC/1.0/" target="_blank" rel="noreferrer noopener">In Copyright - Non-Commercial Use Permitted</a> Format The file format, physical medium, or dimensions of the resource application/pdf Language A language of the resource English Is Part Of A related resource in which the described resource is physically or logically included. Global Ecology and Biogeography Type The nature or genre of the resource Article