# Integrating microclimatic variation in phenological responses to climate change: a 28-year study in a hibernating mammal. Running head: Phenology and microclimate

Anouch Tamian, Vincent A Viblanc, F. Stephen Dobson, Peter Neuhaus, Tracey L Hammer, Anna P Nesterova, Shirley Raveh, Amy L Skibiel, David Broussard, Theodore G Manno, et al.

► To cite this version:

Anouch Tamian, Vincent A Viblanc, F. Stephen Dobson, Peter Neuhaus, Tracey L Hammer, et al. Integrating microclimatic variation in phenological responses to climate change: a 28-year study in a hibernating mammal. Running head: Phenology and microclimate. Ecosphere, 2022, 13 (5), pp.e4059. 10.1002/ecs2.4059. hal-03549269

HAL Id: hal-03549269

https://hal.science/hal-03549269

Submitted on 31 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Integrating microclimatic variation in phenological responses to climate change: a 28-year study in a hibernating mammal.

Running head: Phenology and microclimate

Anouch TAMIAN¹, Vincent A VIBLANC¹*, F Stephen DOBSON¹,², Peter NEUHAUS³, Tracey L HAMMER¹,³, Anna P NESTEROVA⁴, Shirley RAVEH⁵, Amy L SKIBIEL⁶, David BROUSSARD⁷, Theodore G MANNO⁸, Nandini RAJAMANI⁹, Claire SARAUX¹*

1 ¹Institut Pluridisciplinaire Hubert Curien, CNRS, Département Ecologie, Physiologie et

2 Ethologie, 23 Rue du Loess, 67037 Strasbourg, France

3 ²Department of Biological Sciences, Auburn University, Auburn, AL 36849, USA

4 ³Department of Biological Sciences, University of Calgary, Canada

5 ⁴FaunaStats, 16 avenue de l'Europe, Immeuble SXB1, 67300 Schiltigheim, France

6 ⁵Institute of Biodiversity, Animal Health and Comparative Medicine, University of Glasgow,
Glasgow, UK

8 ⁶Department of Animal, Veterinary and Food Sciences, University of Idaho, Moscow, USA

9 ⁷Department of Biology, Lycoming College, Williamsport, PA 17701, USA

10 ⁸Science Department, Catalina Foothills High School, 4300 E. Sunrise Drive, Tucson, AZ

11 85718 USA

12 ⁹Indian Institute of Science Education and Research, Tirupati Andhra Pradesh, India

13

14 * Co-senior authors
15 Author contribution: Conception and design: AT, CS, VAV, PN & FSD; data acquisition:

16 FSD, PN, DB, TLH, TGM, APN, SR, ALS, NR, VAV, CS; data analyses: AT; manuscript

17 preparation: AT. All authors read and commented on the final manuscript.

18

19 Corresponding author: Anouch Tamian, Université de Strasbourg, Centre National de la

20 Recherche Scientifique, Institut Pluridisciplinaire Hubert Curien, Unité Mixte de Recherche

21 7178, 67087 Strasbourg, France. anouch.tamian@iphc.cnrs.fr

22 Open Research statement: Data supporting this publication will be made permanently

23 available on Figshare or Mendeley data upon acceptance of the paper.
24 ABSTRACT

Phenological shifts associated with directional changes in climate, resulting in earlier spring activities, have been documented in several animal species. However, the extent to which species respond to overall climate change versus local climate variation is rarely studied. In addition, climate data are usually averaged over large spatial scales, even though local heterogeneity in habitats may be high, and species might be more susceptible to changes in local rather than global climate conditions. In this study, we examined the effects of spatio-temporal climate variation and climate change on the phenology of a hibernating mountain rodent, the Columbian ground squirrel (*Urocitellus columbianus*). Over 28 years of research (1992-2019), we studied the relationship between the microclimatic conditions experienced by adult and juvenile ground squirrels from four neighboring meadows, and their dates of emergence from hibernation. We used a microclimate model to calculate microclimate variables (local snow depth, soil temperature, air temperature, wind speed, and humidity) at an hourly scale, a 5-m spatial resolution, and at animal height on the study sites over 28 years. Emergence dates varied with age and sex, among years, as well as among and within meadows, with some areas averaging earlier emergence from hibernation than others (up to 10 days of difference). While emergence dates tended to be delayed throughout the study period, long-term temporal changes and interannual variability in emergence dates differed among meadows and depended on individual age and sex. Dates of hibernation emergence were correlated with local climate variables considered either during hibernation or during the preceding summer. Ground squirrels emerged earlier in years or at locations when/where snow melted earlier (years: all individuals excluding 2-year-old males, locations: yearlings and older females), and when the previous summer was less windy (≥3-year-old individuals) and more humid (2-year-old males). Two-year-old male ground squirrels also emerged later in locations where snow depth during winter was higher. Using a microclimate model
49 provided more realistic predictions of phenological responses to climate, highlighting its

50 potential for research on animal responses to abiotic change.

51 **Keywords**: Phenology, climate variability, hibernation, exogenous effects, microclimate

52 model, spatial variability, snow cover.
# INTRODUCTION

The extent to which animals are able to adapt to current climate change (e.g. earlier springs in the temperate zones) and increasing fluctuations in extreme events (e.g. storms, droughts, late snow falls) is of major concern in terms of biodiversity loss and ecosystem function (Walther et al. 2002, Martay et al. 2017, Pires et al. 2018). In particular, variations in winter and spring onsets are occurring globally in seasonal ecosystems around the world (Menzel and Fabian 1999, Kunkel et al. 2004, Kreyling and Henry 2011, Pachauri et al. 2014a, 2014b), affecting the annual timing, known as phenology, of life cycle events (e.g. reproduction) in many animal species (Charmantier et al., 2008; Cohen et al., 2018; Sparks & Yates, 1997). Phenological shifts are generally consistent with the direction of climate change in most animal species, such as earlier spring activities following earlier spring onsets (Parmesan and Yohe 2003, Radchuk et al. 2019). The strength of phenological shifts, however, varies greatly within and among species, between trophic levels, and between geographical locations, from little or no changes to marked changes in seasonal timing (Radchuk et al., 2019; Strode, 2003; Visser & Holleman, 2001; Winder & Schindler, 2004; Chmura et al., 2019; Parmesan, 2007; Thackeray et al., 2016; see Parmesan & Yohe, 2003; Root et al., 2003 for a meta-analyses).

Because global climate change is altering phenological cycles (Menzel 2002, Cleland et al. 2007, Chmura et al. 2019), and because phenological shifts have marked consequences for individual reproduction (e.g. match-mismatch of reproduction and peak resource availability; Durant et al., 2007; Visser & Both, 2005) and population dynamics (Ozgul et al., 2010), it is critical to understand how species phenology responds to fluctuations in climate. The need to distinguish phenological responses to long-term climate change from responses to short-term climate variability complicates this understanding. While climate variability usually considers year-to-year variation in climate, climate change only refers to those
variations that persist for longer periods, usually greater than a decade (World Meteorological Organization, 2015). Another layer of complexity is added by the temporal and spatial scales over which inferences are drawn. Numerous studies rely on satellite or weather station data, often collected several kilometers or tens of kilometers from studied populations (Graae et al. 2012, Potter et al. 2013). However, individuals are often more affected by local climate conditions, which are shaped by spatial heterogeneity in habitat topography, soil composition, or sun exposure (Suggitt et al. 2011, Maclean et al. 2017, Zellweger et al. 2019). This is especially true in rugged areas, such as mountain habitats, where heterogeneity is high over relatively small spatial scales (Huggett 1995, Nagy and Grabherr 2009, Gultepe 2015). Local contrasts in slope angles, sun exposure, soil composition, and vegetation cover can cause climate to vary, through variations in solar radiation, wind exposure or humidity, from scales of 1-10 kilometers (i.e., mesoclimates) to scales of 0-100 meters (i.e., microclimates) (Barry 1992). Due to a lack of climate data at a fine resolution (<1 km), few studies have examined phenological responses to microclimatic conditions (James et al. 1984, Suggitt et al. 2011, Potter et al. 2013, Varner and Dearing 2014, Paruchuri et al. 2019).

In this study, we examined the relationship between the phenology of hibernation and climate in a mountain hibernator, the Columbian ground squirrel (*Urocitellus columbianus*). We focused on hibernation emergence dates (ED) for individual ground squirrels from four neighboring meadows monitored over a period ranging from 18-28 years, depending on the meadow. Our objectives were to (1) examine spatial (among and within the four meadows) and temporal (over the years) variation of ground squirrel ED and to (2) determine how this variation might be explained by local and temporal variation in microclimatic conditions (scale of tens of meters, as experienced by the animals). Columbian ground squirrels are colonial ground-dwelling sciurid rodents inhabiting burrow systems in subalpine meadows of the northern regions of the Rocky Mountains. Columbian ground squirrels are an interesting
model system to investigate changes in phenology in relation to microclimates. These relatively long-lived animals (the oldest squirrel in our data set lived up to 14-yo.) hibermate for ~8 months of the year, and are only active for a short summer period during which reproduction occurs (Murie and Harris 1982, Young 1990b, Dobson et al. 1992). Because of this short active season, phenological shifts may be important in determining animal reproduction or ability to fatten for winter. A long-term delay adult female ED, correlating with lower spring temperature and delayed snowmelt, has already been reported (9.4 days over a 20-year period) in one of our four populations (Lane et al., 2012). While this study suggested phenological plasticity in ground squirrel ED, likely due to climate change (see also Dobson et al., 2016), it only considered adult females from a single population and used wide-scale climate variables to test for their effects on ground squirrel ED (the weather data came from a meteorological station located in a prairie habitat, some 55km from the study site).

Yet, although some emigration occurs (Wiggett et al., 1989), adult ground squirrels are largely highly philopatric and have a relatively small ranges over their lifetimes (e.g. female range rarely exceeds 1 km² and male range around 4.2km²; Arnaud et al., 2012; Festa-Bianchet & Boag, 1982; Harris & Murie, 1984; Murie & Harris, 1978). Thus, these animals are more likely to respond to very local climate conditions than to wide-scale climate variations. In addition, previous studies have shown that adult males emerge from hibernation first, followed by adult females, and lastly yearlings (Murie and Harris 1982, Dobson et al. 1992). However, data was lacking on whether different age or sex-categories show similar or different phenological plasticity to local climate, both within and between neighboring populations (resolution of ~300m).

To study the relationship between emergence date and climate, we performed a three-step analysis. First, we analyzed how variable ED was both within and among four different meadows over a 28-year period. We analyzed variation in ED in relation to individual (sex
and age) and environmental characteristics (year, hibernation location) that were expected to affect hibernation emergence. Second, although our four study sites are located in the same valley, they present observable variability in slope, topography, and exposure to sun both within and among sites. Thus, we expected local climate variations to occur, affecting variation in ED both among meadows and within different parts of a given meadow. In the absence of local weather stations, we used a microclimate model (see Kearney, 2020; Kearney et al., 2020; Lembrechts & Lenoir, 2020; Maclean, 2020; Saleeba et al., 2020) to predict climate since 1992 on the four study sites, at a 5-m resolution, and at the height of the animal (10 cm above ground, or 1-m below ground, roughly the depth of a hibernaculum; Young, 1990). We then analyzed climate variation and tested whether we could detect a directional change in local climate patterns over the past 28 years (viz. climate change). Finally, we tested whether variation in ED could be explained by the local climate conditions encountered by the animals. Because emergence from hibernation is strongly related to the energetic state of individuals at the end of the hibernation (Norquay and Willis 2014, Williams et al. 2014, Fietz et al. 2020), there are two theoretical ways that climate could affect ground squirrel emergence dates: 1) by modifying individual energy expenditure during hibernation (Humphries et al. 2002, Pretzlaff and Dausmann 2012, Geiser 2013), and 2) through carry-over effects from one active season to another (Norris 2005, Saino et al. 2017, Moore and Martin 2019). On one hand, the energy available and stored by individuals during the active period should affect their ability to fast during hibernation (as Columbian ground squirrels are "fat-storing" hibernators, Murie and Boag 1984, Humphries et al. 2003). On the other hand, harsh climate conditions during winter could also affect ground squirrel energy expenditure, and thus influence individual ability to fast and hibernate (Davis 1976, Humphries et al. 2004, Fietz et al. 2020).
As the environment is multifactorial by nature, our analyses included several climate variables that could affect either energy storage before hibernation or energy expenditure during hibernation (as recommended by McGuire et al. 2021). Specifically, we expected areas or years with earlier snowmelt and higher underground temperature to be associated with earlier ED (see Lane et al. 2012). The presence and depth of the snow cover directly impacts soil insulation and ground temperature, potential external cues for hibernating individuals. In addition, we expected carry-over effects of climate on ED to occur from one year to the next. During hibernation, ground squirrels rely almost exclusively on body fat (Shaw 1926).

Overwinter survival depends on fat reserves accumulated during the previous season (e.g. Murie & Boag, 1984), and low fat deposition is likely to affect ED the subsequent year. We did not measure body mass before hibernation in our study, but we reasonably expected foraging conditions (length and quality of forage) during the previous summer to directly reflect the amount of fat accumulated and burned during winter (Dobson et al. 1992, Neuhaus 2000). Thus, we expected unfavorable forage conditions – e.g. low rainfall (reflected by humidity) to be associated with low primary productivity (Dobson and Kjelgaard 1985, Neuhaus et al. 1999, Dobson et al. 2016), and strong winds to be associated with more time spent by animals in vigilance rather than foraging (Fairbanks & Dobson, 2007). Through the decrease of fur thermal insulation (McCafferty et al. 2017), humidity might also negatively affect the time foraging, which will have detrimental effects on the accumulation of fat stores and ultimately on hibernation duration. This may cause individuals to emerge earlier due to the early depletion of overwintering fat stores, or conversely to delay emergence by more pronounced hypometabolism and a longer torpor period until resource availability is high again (i.e. avoiding early arousals when resource availability is low). Juveniles, which are less able to mitigate the survival costs of low-fat energy stores during hibernation, might experience this more than adults (Murie & Boag, 1984). Alternately, more favorable
conditions during the previous summer could also allow for longer hibernation periods and later ED in juveniles, allowing them to avoid the period of high territorial aggression displayed by reproducing adults early in the season (Murie & Harris, 1978; Murie & Harris, 2011). Taken together, the joint investigation of summer and winter condition effects on ED in different neighboring populations was expected to provide us with a comprehensive understanding of how these mountain-adapted rodents deal with local heterogeneity in weather conditions, and whether different populations respond differently to long-term changes in climate patterns.

## MATERIALS AND METHODS

### Emergence date sampling

#### Study site and population monitoring

Columbian ground squirrels were monitored in the Sheep River Provincial Park in Alberta, Canada (50°38′N, 114°39′W), located in the foothills of the Rocky Mountains (Fig. 1.A). This study area includes four different meadows surrounded by mixed forests (primarily lodgepole pine, white spruce, quacking aspen, and birch trees) and composed of vegetation that is representative of grasslands in montane sub-regions (Alberta Parks 2008). We monitored ground squirrel colonies on: Meadow A (3.6 ha, from 1510-1535m, monitored since 1994), Meadow B (2.6 ha, from 1500-1540m, monitored since 1992, meadow on which the Lane et al. 2012 study was done), Meadow C (1.7 ha, from 1540-1560m, monitored since 1999) and Meadow DOT (5.7 ha, from 1527-1570m, monitored since 2001). Meadows A, B, and C slopes are southeast facing, whereas Meadow DOT slope is southwest facing.

We collected individual data with comparable methods on all four meadows.

Individuals were trapped and permanently marked as pups, or when they first appeared in the
meadow as immigrant adults, with application of unique metal numbered ear tags (Model no. 199
1, National Band & Tag Co.). Sex was determined based on visual inspection of genitalia 200
(Murie & Harris, 1982). Each year, the meadows were monitored daily from before the first 201
emergences from hibernation (early April each year) to the end of lactation for breeding 202
females. Unfortunately, dates of entry into hibernation were not recorded and we were not 203
able to follow the entire cycle of annual phenology. From mid-April to the end of May each 204
year, ground squirrels were captured within a couple of days of emergence from hibernation 205
(usually on the day of emergence or the following morning) using live-traps (National Live 206
Traps; Tomahawk Co., WI, USA: 13 x 13 x 40 cm³) baited with peanut butter. The first day of 207
observation or capture (usually the same) of a given squirrel was used as a proxy for 208
hibernation emergence date (ED), typically confirmed by the squirrel's appearance and 209
physical condition (large skin flakes from hibernation lost within few days, large abdominal 210
flaps of skin where fat reserves were accumulated and lost, no defecation upon capture; Murie 211
& Harris, 1982). Emergence was recorded as an ordinal date (i.e. number of days after 1st 212
January each year). Emergence data were not collected in 2006 on Meadow DOT, in 2003- 213
2006 on MA, and in 2002-2003 on MC. The location of each initial capture was recorded at a 214
± 1.0m resolution, on 10 x 10 m Cartesian grids of color-flagged wires at each meadow. 215
Hibernacula are usually used by the animals during a few days after emergence until 216
collective burrows were re-excavated for the active season. First capture locations should thus 217
provide a reasonable approximation of the area where individuals hibernated. 218

Recorded Cartesian coordinates were converted to GPS coordinates using a database 219
of GPS and Cartesian coordinates collected in 2018-2019. We estimated the error associated 220
with the conversion of Cartesian to GPS coordinates at a median of 5.5 ± (SD) 3.4m (n = 225 221
coordinates). Coordinates were collected for 74% of the emergences, 4351/5873. The missing 222
emergence sites were before 1994 on Meadow A, 1999 on Meadow B, 2006 on Meadow C, 223
and 2012 on Meadow DOT. Emergence locations were also not recorded in 2000-2007 on MA. Thus, sample sizes varied among years.

Age classification

Because male Columbian ground squirrels emigrate, and new individuals occasionally leave and arrive on meadows, some individuals on the Meadows A, B, C, and DOT were not known from birth (19% of the cases overall, 455/2436). For these individuals, recorded age resulted from estimates of physical characteristics mainly distinguishing between yearlings, 2-year-old, and older individuals. To establish meaningful age bins in further analyses we conducted a preliminary analysis investigating how emergence date varied according to age using individuals of known age (n = 4681 observations, N = 1981 individuals). We ran a linear mixed effects model (LMM) of emergence date (dependent variable) as a function of age (independent categorical variable), specifying individual identity as a random factor in the model to account for repeated measures of individuals among years. Tukey HSD tests were applied for post-hoc comparisons. Emergence date of one-year-old individuals (ordinal date ± SE = 126.5 ± 0.3 days) was significantly later than all other age categories (Tukey HSD, p < 0.001). For 2-year-old individuals, emergence ordinal date (121.7 ± 0.3 days) was significantly earlier than yearlings (p < 0.001) but significantly later than older age categories (116.7 ± 0.2 days; p < 0.001). No other significant differences were found among older age categories (all p > 0.05). Thus, individuals aged 3 and older were combined in a unique age category (≥3-year-old) for subsequent analyses allowing us to use all existing ED data (n = 5873 observations in total: N = 2436 individuals). Sample sizes per meadow, year and age and sex classes are provided in Appendix S1 (Table S1).
247 To evaluate spatial heterogeneity in emergence dates within a meadow, we
248 investigated the spatial autocorrelation of emergence dates. Specifically, we examined
249 whether correlation between emergence dates at two locations depended on the distance
250 between them. To do so, we used variograms (Matheron 1963) that measured the value of the
251 spatial variance at different distance intervals. For any two data points, we measured their
252 spatial distance and estimated the variance in ED. The variance was then averaged over bins
253 of distance intervals. Empirical variograms were based on all emergence locations, but with
254 pairs of points considered only if points were in the same year and within the same meadow.
255 We modelled variograms using a least square method (Rivoirard et al. 2008). The empirical
256 variogram of emergence date exhibited spatial structure (i.e., increasing variance when the
257 distance between sampling points increased before levelling off; see Appendix S1, Figure
258 S1), meaning that the emergence dates were not random, such that spatially close locations
259 had closer emergence dates than more distant locations. Based on the modelled variogram, we
260 found that variance in emergence date stabilized at a range of about 40 m. Thus, hibernation
261 emergence locations were grouped into spatially relevant clusters according to their GPS
262 coordinates to understand whether individuals on some parts of the meadow emerged earlier
263 than others.

264 We built clusters by grouping emergence locations according to an iterative algorithm
265 (adapted from Woillez et al., 2007, see also Saraux et al., 2014). The algorithm started with
266 each location as its own cluster. Each iteration of the algorithm then: (1) calculated the
267 distances between the center of gravity of all clusters on the meadow; (2) merged the two
268 closest clusters if the shortest distance among all cluster distances was below a threshold
269 distance of 40 m (according to above variogram results). The process was reiterated until
270 none of the distances between clusters were lower than 40 m. In total, we defined eleven,
eight, seven, and twelve clusters for Meadows A, B, C, and DOT respectively (named differently depending on the meadow), which corresponded to hotspots of ground squirrel emergences in the habitat (compared to zones where ground squirrels were rarely seen emerging). Emergence locations and built clusters were imported into QGIS (3.10.12 long-term released version, QGIS.org, QGIS Geographic Information System, 2021) and presented on a Google satellite map © 2021 (see Fig. 1.B).

## 277 Microclimate data

We calculated hourly microclimatic conditions using the Kearney et al. (2020) microclimate model, at a fine-scale resolution (5 x 5 m grid) for each meadow during the entire monitoring period (from 1992-2019). The modelling function (*micro_ncep* from NicheMapR package, Kearney et al., 2020) downscales global atmospheric climate-forcing data (NCEP Reanalysis, historical data from the National Centre for Environmental Predictions, Kalnay et al., 1996, provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, from their Web site at https://psl.noaa.gov/), using terrain-specific corrections, including spatial variations in slopes (from the Mapzen elevation and Mapzen terrain service, via 'elevatr' package in R; Hollister et al., 2020), soil composition and hydrological properties (acquired from the 'SoilGrids' database in R; Hengl et al., 2017), orientation, and hill shade. The microclimate model increases the quality and resolution of large-scale weather data, and has been empirically validated (Kearney, 2020; Lembrechts & Lenoir, 2020; Maclean, 2020; Saleeba et al., 2020). Because Columbian ground squirrels live in open grass meadows with adjacent trees (see Fig. 1), we estimated climatic variables without specifying shade due to vegetation cover.

Our model estimated snow depths (cm) on the ground as well as temperatures (°C), wind speed (m/s), and relative humidity (%) at 10 cm above ground (roughly the height of a
ground squirrel). In addition, we estimated temperature at various soil depths below ground (0, 2.5, 5, 25, 50, 75, 100, 125, 150, and 200 cm). Hibernating burrows are usually around 70 cm in depth and rarely exceed 1 m (Young, 1990). Soil temperatures at various depths were highly correlated (see Appendix S1, Figure S2). Thus, we used 1-meter-deep temperatures as representative of hibernaculum conditions.

There was no reasonably close weather station with records of climate date for the entire study period (the closest one being 58 km away at lower altitude in the plains). In addition, weather stations record temperature and precipitation variables, at a standard height of 2 m, which is not representative of what a ground-squirrel might experience a few centimeters above ground. Thus, using the microclimate model enabled us to obtain climate data over the entire study period at a high resolution and at the height of the animal, and further allowed us to measure snow cover (see Appendix S1, Figure S3 for a validation of the microclimate model).

To estimate spatial and temporal variation of climate, we averaged the model-generated time series at two different spatial scales for each meadow: i) over the entire meadow and ii) separately within each cluster. Hourly data were averaged over 24-hour periods. Based on the average ordinal day of first hibernation emergence (older males) across all years (15 April ± 2 days) and previous studies showing that immergence occurs close to the beginning of August (Young 1990b, Neuhaus 2000), we further divided the year in two phases: the active season (15 April – 1 August), and the inactive (hibernation) season for the remainder of the year. This allowed us to test for possible effects of climate conditions during the winter and the previous summer on individual emergence dates. Besides air temperature (°C), soil temperature at 1m depth (°C), relative humidity (%), wind speed (m/s), and snow depth (cm), we considered two additional climate variables: first day of the calendar year without snow cover on the ground (over the entire meadow or within each given spatial
cluster, in ordinal date), and number of days with snow falls during the active period for
investigation of late snowstorms on ED.

322 **Data analyses**

323 *Quantifying variance in emergence date*

324 *Emergence date variance partitioning*

325 To identify how much variance in ground squirrel emergence dates (ordinal days) was
326 explained by endogenous (sex, age category, individual ID), spatial (meadow, emergence
327 cluster within the meadow), and temporal (year) factors, we used a linear mixed effects model
328 (LMM) with all variables specified as random terms (with 1|variable):

329 *Emergence date* ~ (1|sex) + (1|age<sub>category</sub>) + (1|ID) + (1|meadow/cluster) + (1|year) [1]

330 This model allowed us to partition variance components between the different random terms;
331 (1|meadow/cluster) representing both the effect of the meadow, and the effect of the cluster
332 nested within the meadow.

333 *Effects of sex and age on emergence date*

334 Because sex and different age (especially sexually mature versus immature
335 individuals) classes face different reproductive constraints, we tested the effects of individual
336 endogenous characteristics (sex and age specified as independent variables) on emergence
337 date, while controlling for spatial and temporal effects (set as random terms):

338 *Emergence date* ~ sex + age<sub>category</sub> + sex*age<sub>category</sub> + (1|ID) + (1|meadow/cluster) + (1|year)

339 [2]
340 *Spatio-temporal variation in emergence date*

341 *Spatial variation in emergence date*

342 Because of differences in slope and exposure at different meadows, we expected
343 emergence dates to vary spatially. To test for spatial variation in emergence date, we assessed
344 among and within meadow differences using separate LMMs. Differences were tested by
345 specifying emergence date as the dependent variable, meadow (or cluster) as an independent
346 variable, and year and the sex-age category (determined based on sex and age results) as
347 random terms:

348 *Emergence date ~ meadow*[MA vs. MB vs. MC vs. DOT] + (1|year) + (1|sex-age<sub>category</sub>) + (1|ID) [3]*

349 *Emergence date ~ cluster* + (1|year) + (1|sex-age<sub>category</sub>) + (1|ID) [4] (one LMM per meadow)
350

351 *Temporal variation in emergence date*

352 Individuals might be phenotypically plastic, and time their ED to interannual
353 variations, including directional changes in climate variables. We differentiate between
354 interannual variation in ED and long-term directional changes in ED. Because emergence
355 locations were not always known (especially in the early years, see above), we did not include
356 spatial clusters in the following models [5] and [6], enabling us to study potential changes in
357 time over a larger data set (e.g. starting in 1992 instead of 1999).

358 **Interannual variation in EDs**

359 To explore interannual variation in hibernation emergence dates, we built generalized
360 additive mixed models (GAMMs) with emergence date as the response variable, and year as
361 the independent variable. This facilitated identification of changes between years, and
362 detection of nonlinear temporal variation. Because temporal trends could differ among
363 meadows or according to the individual sex-age category, we considered interactions
364 (specified with “by =”) between the continuous variable (year) and these factors by
365 constructing separate GAMMs, keeping individual ID as a random variable:

366 *Emergence date* ~ *s(year, by = meadow)* + *s(year, by = sex-age<sub>category</sub>)* + *sex-age<sub>category</sub>* +
367 *meadow* + (1|*ID*) [5]

368 **Long-term directional changes in ED**

369 In a context of climate change, most phenological studies showed advanced spring
370 season (Walther et al. 2002, Cleland et al. 2007, Dingemanse and Kalkman 2008). Yet, in
371 Columbian ground squirrels, females older than 3 displayed a delay in emergence over a 20-
372 year period (Lane et al. 2012). To test whether this long-term directional pattern was
373 consistent over the 4 populations and over other sex/age categories, we ran the previous
374 model (model [5]) as a LMM, and assessed interactions between continuous (year) and
375 categorical (meadow or sex-age<sub>category</sub>) variables:

376 *Emergence date* ~ *year* + *sex-age<sub>category</sub>* + *meadow* + *year*\**meadow* + *year*\**sex-age<sub>category</sub>* +
377 (1|*ID*) [6]

378 *Spatio-temporal variation in climate*

379 *Spatial variation in climate variables*

380 To investigate microclimates, we assessed correlations among the 11 previously
381 defined climate variables (air temperature during the active and inactive period, soil
382 temperature at 1m depth during the active and inactive period, relative humidity during the
383 active and inactive period, wind speed during the active and inactive period, snow depth
384 during inactive period, first day of calendar year without snow and number of days with snow
385 during the active period), using a Principal Component Analysis (PCA) on data averaged by
386 year and by cluster on each meadow. The objective was to identify principal variations related
387 to spatial variables (clusters and meadows) given the large amount of climate data at a micro-
scale within the meadows. Spatial differences in climate were analyzed from resulting

components (PCs), through LMs at two different scales:

PCs ~ meadow [7]

PCs ~ cluster [8]

Using our results (see below), only climate variables expected to affect emergence

date were kept for subsequent analyses. These included (1) snow depth (depth<sub>snow</sub>) and soil

temperature (temp<sub>soil</sub>) at a 1-meter depth during hibernation (inactive phase), and (2) air

temperature (temp<sub>air</sub>), relative humidity and wind speed during the previous summer (active

phase), in addition to the ordinal date of spring snowmelt (day<sub>snowmelt</sub>) and the number of days

with snow during the active period (days<sub>snow</sub>).

Temporal variation in climate variables

Because interannual variability in climate often result from different processes: trend

in climate (e.g., due to climate change), cyclic effects (e.g., El-Niño) and random year-to-year

variation, we decided to use a 2-step approach. First, we focused on short-term variation using

GAMs. Then we looked for directional trend in ED (as expected as a result of climate change)

using LMs.

Climate variation

To explore how climate varied, we ran GAMs for each climate variable as a function

of year. This allowed identification of changes between years, and detection of nonlinear

temporal variation. Because the interannual variability of climate could differ among

meadows (or among clusters), we also considered interactions of years with these factors in

the following GAMs:

Climate variable ~ s(year, by = meadow) + meadow [9]

Climate variable ~ s(year, by = cluster) + cluster [10] (one GAM per meadow)
412 Climate change

413 To test for long-term directional climate changes over 28 years, we ran the previous

414 GAMs [9] & [10] as LMs [11] & [12], with each climate variable as a dependent variable and
415 the year, meadow (or cluster) and their interaction as independent variables:

416 Climate variable ~ year*meadow + meadow [11]

417 Climate variable ~ year*cluster + cluster [12] (one LM per meadow)

418 Climate – emergence date relationships

419 We tested how climate might affect emergence date through two processes: energy
420 expenditure and carry-over effects of energy stored before hibernation. Because Columbian
421 ground squirrels start breeding at 2-years-old, the emergence date of ≥ 2-year-old individuals
422 may vary with climate differently between males and females. Thus, we considered five sex-
423 age categories for these analyses: yearlings, 2-year-old females, 2-year-old males, ≥3-year-
424 old females, and ≥3-year-old males.

425 Mean relationships between climate and emergence date

426 To test for overall relationships between microclimates and hibernation emergence
427 dates, emergence dates were averaged by year, meadow, and sex-age category. We built
428 separate LMMs for the different age-sex categories, with emergence date as the dependent
429 variable and climate variables (averaged by meadow and year) as independent variables:

430 Emergence date ~ depthsnow_winter + tempsoil_winter + tempair_active n-1 +
431 humidityactive n-1 + windactive n-1 + daysnowmelt + dayssnow_active + (1|Year) +
432 (1|Meadow) [13] (one LMM per age-sex category)

433 Individual responses in emergence date might differ due to different constraints
434 between spatial heterogeneity (related to spatial living locations) and temporal heterogeneity
435 (plasticity that matches interannual variation). To distinguish between temporal vs. spatial
effects of climate on hibernation emergence dates, we ran two further analyses, building separate models for different sex-age categories.

## Year-centered approach: are spatial differences in ED related to spatial variation in climate?

To investigate spatial effects of climate variables on emergence date, we removed all temporal variation by centering emergence date and climate variables by year and sex-age category (mean of the years and category subtracted to each observation). During the active season, ground squirrels might visit most locations on the meadow in a single day. Because spatial clusters were defined from hibernation emergence locations and not the total active range of individuals, we included in the model only climate variables considered during hibernation (inactive phase). For each individual emergence date, we associated the climatic conditions of the given winter averaged over the spatial cluster in which the individual emerged. Thus, the linear mixed model was:

$$ \text{Emergence date} \sim \text{depth}_{\text{snow\_winter}} + \text{temp}_{\text{soil\_winter}} + \text{day}_{\text{snowmelt}} + \text{days}_{\text{snow\_active}} + (1|\text{cluster}) + (1|\text{ID}) \text{ [14]} $$

## Spatial-centered approach: are interannual differences in ED related to interannual variation in climate?

To investigate temporal effects of climate variables on emergence date, we centered emergence date and all climate variables (including those considered during the active phase) by cluster and sex-age category. This way, spatial variations in climate were entirely removed from the data, and we tested whether yearly differences in climate could explain differences in emergence date among years.

$$ \text{Emergence date} \sim \text{depth}_{\text{snow\_winter}} + \text{temp}_{\text{soil\_winter}} + \text{temp}_{\text{air\_activen-1}} + \text{humidity}_{\text{activen-1}} + \text{wind}_{\text{activen-1}} + \text{day}_{\text{snowmelt}} + \text{days}_{\text{snow\_active}} + (1|\text{year}) + (1|\text{ID}) \text{ [15]} $$
Statistical analyses were done in R v 4.0.3 (R Core Team, 2019). All tests were conducted with a probability error threshold of 5%. Results presented are means ± SE, along with the number of observations (n) and the number of individuals (N). Where appropriate (i.e. for LMMs and LMs), we ensured model residuals were normally distributed by visual inspection of density distributions, Q-Q plots, cumulative distribution functions, and P-P plots using the “fitdistrplus” package in R (Delignette-Muller and Dutang 2015). When the assumptions of a normal distribution of residuals were not met, we applied a Box-Cox transformation (from MASS R package, Box & Cox, 1964) to the data (monotonous transformation), or removed extreme values from the data (i.e., emergence dates later than the 22nd of May for which we likely missed the actual emergence date, representing 2.7% (116/4351) of observations). Results with and without transformation or extreme values yielded similar results, and we present results obtained from complete and untransformed data sets for clarity (exceptions are specified in the results). Post-hoc multiple comparisons were run to compare all factor levels of categorical variables, using the Least Square Means (LSM) test, to account for multiple testing (Lenth 2016). GAM(M)s were run using the ‘gam’ and ‘gamm’ function in R (‘mgcv’ package, Wood, 2017). Principal Component Analysis was run using the ‘PCA’ function, and components accounting for a meaningful amount of variance (i.e. from 75%, O’Rourke & Hatcher, 2013) were then considered in the models [7] and [8] (‘FactoMineR’ package, Lê et al., 2008). Whenever interaction effects between independent variables were included in LMs and LMMs, we compared the models with and without the interactions, and chose the final model according to Akaike’s Information Criterion (AIC) (lowest AIC and for ΔAIC < 2, we chose the most parsimonious model, Burnham & Anderson, 2004). However, we also analyzed meadows separately to assess difference in trends between populations. The results on the link between climate and emergence date are
presented from the full models (i.e. with all tested variables included). Multicollinearity was checked and when variables were collinear, the one with the highest variance inflation factor (VIF) was removed, meaning that only variables with variance inflation factors below three remained in the models (Zuur et al. 2010). Finally, to allow comparisons of effect sizes in the last models [13], [14], and [15], we standardized (by subtracting the mean and dividing by the standard deviation) all explanatory climate factors.

## 490 RESULTS

Quantifying variance in emergence date

Emergence date variance partitioning

Across all individuals, meadows, and years, mean emergence date from hibernation was 2 May ± 4 days (ordinal date 122 ± 4) and ranged from 6 April to the 30 May (ordinal date 96-150). The model [1] with age, sex, meadow, cluster within meadow, year, and individual as random effects explained 65.7% of the total variance in emergence date (LMM, n = 4351 data points). Variables contributing most were age (26.3%), year (16.1%), and individual identity (14.6%). Within-meadow cluster, sex, and meadow explained fairly trivial amounts of variation, 3.9%, 3.1%, and 1.7% of the variance in ED, respectively.

Effects of sex and age on emergence date

Age, sex, and their interaction, were all retained in the best model explaining emergence date (LMM [2], n = 4351, N = 1950 individuals, 26 years). We found differences in emergence date between males and females only in individuals ≥3 years old, with males emerging around 7.3 ± 0.4 days earlier than females (post-hoc LSM, z = -18.3, p < 0.001). Post-hoc comparisons revealed significant differences in emergence dates between four categories: yearlings, 2-year-old individuals, ≥3-year-old females and ≥3-year-old males (LSM, all p < 0.001). On average, older males (≥3-year-old) emerged first around 21 April
508 (111.2 ± 0.2 days), followed by older females (≥3-year-old) around 28 April (118.5 ± 0.2 days), 2-year-old individuals around 2 May (121.8 ± 0.3 days), and finally one-year old individuals around 6 May (125.5 ± 0.2 days).

511 *Spatio-temporal variation in emergence date*

512 *Spatial variation in emergence date*

513 Analyses on spatial differences in hibernation emergence dates showed that
514 individuals on Meadow DOT emerged 2.2 ± 0.5 days and 1.8 ± 0.4 days earlier than
515 individuals on Meadow C and Meadow B, respectively. Meadow C and Meadow B
516 individuals emerged 3.3 ± 0.4 days and 3.6 ± 0.3 days earlier than Meadow A individuals
517 (LMM [3], post-hoc LSM: all p < 0.001, n = 5873, N = 2436 individuals, 28 years). Similarly,
518 within-meadow differences in emergence dates were found between different clusters (see
519 cluster colors on Fig. 1.B). As an example, on Meadow B, individuals in “bottom-west”
520 cluster emerged between 1.9 ± 0.5 and 8.8 ± 2.5 days earlier than individuals from other
521 clusters (LMM [4], post-hoc LSM: p < 0.001, n = 1315, N = 579 individuals, 21 years). In
522 contrast, individuals from the “top” cluster emerged between 4.6 ± 0.9 and 8.1 ± 0.8 days
523 after individuals from other clusters (post-hoc LSM: p < 0.001, n = 1315, N = 579
524 individuals, 21 years) (see Fig. 1.B and Appendix S1, Table S2).

525 *Temporal variation in emergence date*

526 **Interannual variability in ED**

527 Emergence date exhibited marked interannual variations over the course of the study
528 on all meadows and for each age-sex category (GAMMs [5]; 2.09 < edf < 8.73, 3.43 < F <
529 31.61, p < 0.001, Fig. 2), except for older males (≥3-year-old) on Meadow B and C (p =
530 0.158 and p = 0.652, respectively). Cyclic patterns over time were observed with a period
531 varying from 6-10 years, depending on the period or the meadow considered. Older males
(≥3-year-old) showed less variation in emergence date over time on Meadow A than on Meadow DOT. Emergence date of yearlings and older females on Meadow A was especially delayed between 1995-97 as compared to 1999, as well as in 2011 compared to 2014-17. On Meadow B, yearlings' and ≥3-year-old females' emergence date exhibited shorter fluctuations at the start of the monitoring period (1992-2002) than in recent years (2005-2019). Emergence date appeared to be delayed from 2005 to 2009, then advanced from 2009-2016, and delayed again from 2016-2019. On Meadow DOT, emergence date for yearlings and older individuals showed lower variation than on other meadows, with fluctuations of smaller amplitude. Meadow C individuals exhibited the greatest variation during the most recent years, emergence dates being advanced between 2013 and 2015 and delayed between 2016 and 2018. In general, the emergence date of 2-year-old individuals showed slightly lower fluctuations and a more stable pattern. On all meadows, individual emergence dates appeared to be earlier in 2015 and 2016 when compared with other years. Whereas the emergence date of Meadow A and Meadow B individuals was later in 2010 than other years, it was earlier on Meadow C and Meadow DOT.

## 547 Long-term directional changes in ED

When testing for long-term directional linear changes in emergence dates, the interaction between year and meadow was retained in the best model (LMM [6], ΔAIC = 35, n = 5873, N = 2436 individuals). On Meadow DOT, emergence dates of individuals did not significantly change in a directional manner through time for any of the sex-age categories (LMM, interaction year*sex-age not retained, ΔAIC = 5, t = 1.094, p = 0.274, n = 1303 observations, N = 503 individuals) (see Appendix S1, Figure S4). On Meadow C, regardless of sex and age (interaction year*sex-age not retained, ΔAIC = 2), emergence from hibernation was significantly delayed over time, averaging a 4-day later emergence over an 18-year period (LMM, +0.22 ± 0.06 day/year, t = 3.769, p = 0.001, n = 704 observations, N = 385
individuals). On other meadows, temporal trends differed according to the sex-age category ($\Delta$AIC = 17 and 22 for meadow A and B, respectively). On Meadow A, over a 22-year period, emergence was advanced by 2 days for 2-year-olds (LM, -0.11 ± 0.05 day/year, t = -2.290, p = 0.023, n = 399 observations, N = 399 individuals), delayed by 3 days for ≥3-year-old individuals (LMMs, ≥3-year-old Females: +0.14 ± 0.04 day/year, t = 3.254, p = 0.001, n = 767 observations, N = 290 individuals; ≥3-year-old Males: +0.13 ± 0.04 day/year, t = 3.058, p = 0.003, n = 368 observations, N = 155 individuals) and stable in yearlings (LM: p = 0.053, t = -1.936, n = 625 observations, N = 625 individuals). On Meadow B, the hibernation emergence of yearlings, 2-year-old individuals and ≥3-year-old females was delayed by 11, 7 and 5 days, respectively over 28 years (LMs; yearlings: +0.40 ± 0.04 day/year, t = 9.110, p < 0.001, n = 631 observations, N = 631 individuals; 2-year-old: +0.26 ± 0.06 day/year, t = 4.678, p < 0.001, n = 318 observations, N = 318 individuals; LMM: ≥3-year-old Females: +0.17 ± 0.05 day/year, t = 3.455, p < 0.001, n = 531 observations, N = 164 individuals), whereas the emergence date of ≥3-year-old males did not significantly change through time (LMM; p = 0.379, n = 227, N = 83).

## Spatio-temporal variation in climate

### Spatial variation in climate variables

The three most influential components of the PCA explained 78.3% of the variance in climate variables and each generated an eigenvalue higher than one. The first component (PC1), explaining 39.4% of the variance, primarily reflected air and soil temperature during the summer (active period), soil temperature during the winter (inactive period), and negatively correlated to wind speed and humidity during the summer (see Fig. 3 and Appendix S1, Table S3). Winter air temperature and humidity were positively correlated with the second component (PC2, explaining 24.5% of the variance), whereas snow depth during winter negatively loaded on PC2. The third component (PC3, explaining 14.4% of the
582 variance), was mainly related to wind speed during the winter. The three principal
583 components varied depending on the meadow (LMs [7], p < 0.001), and depending on the
584 cluster (LMs [8], p < 0.001, see Appendix S1, Figure S5).

585 As presented on Fig. 3.B and 3.D, clusters on Meadow A were between 9.9% and
586 29.2% windier during the inactive period, between 16.3% and 26.4% windier during active
587 periods, and between 0.6% and 2.7% more humid during the active period than clusters from
588 other meadows. Conversely, soil temperatures of clusters on Meadow B were higher, during
589 both winter (between 1.7% and 4.5%) and summer (between 2.4% and 6.6%), and the air
590 during summer was between 1.4% and 4.3% warmer than on other meadows. Despite those
591 general patterns, within-meadow variability was also quite high among clusters, especially on
592 Meadow B, which displayed both very positive and negative PCs clusters (Fig. 3 B & D and
593 Appendix S1, Figure S5).

594 *Temporal variation in climate variables*

595 Climate variation

596 As predicted, all climate variables exhibited temporal variability (GAMs [9] and [10];
597 6.51 < edf < 8.78, 2.79< F < 15.33, all p < 0.01, Fig. 4). Interannual variation in these
598 variables was similar among meadows (no year*meadow interactions retained, N = 28 years),
599 although the absolute values in wind speed and relative humidity varied across meadows. This
600 confirmed previous results that Meadow A was the windiest and most humid meadow, while
601 Meadow B was the least humid and windy, Meadows C and DOT being intermediate.
602 Summer air temperature exhibited particularly high interannual variation, increasing and
603 decreasing from 1 to 4°C from one year to the next (e. g. between 1997 and 1998), fluctuating
604 on average between 7-13°C. During 2006-2007, air temperature was particularly high, and
605 much warmer than the next four years. Conversely, the relative humidity during summer was
606 particularly low in 2006-2007 compared to the next four years, staying relatively stable over
the rest of the study period. Wind speed during the active period showed regular cycles oscillating between 0.8-1.2 m/s every five years. During the inactive period, soil temperature varied between 2-5 °C within approximately five-year periods. Winter snow depth showed interannual variation as well, with especially low values in 2018 and 2019. Similarly, the ordinal day of snowmelt varied between mid-February and 10 April, depending on the year.

Number of days with snow during the active period contrasted between years of no late snowstorms (1992-1995, 1997-1998, 2004-2006, 2012, 2015-2016) and years of late snow falls (between 2000-2003, and between 2008-2013).

## Climate change

Temporal trends in climate did not differ between meadows (LMs: no year*meadow interactions retained, N = 28 years, see Appendix S1, Figure S6). Over 28 years, air temperature increased by 1.40°C during the active season (LM [11], Active: 0.05 ± 0.01°C/year, t = 3.510, p < 0.001). Winter (inactive period) soil temperature at a 1-m depth increased at a slower rate of 0.03 ± 0.01°C/year (i.e. 0.84°C over 28 years, LM, t = 5.642, p < 0.001). Conversely, the snow depth during the inactive period (winter) decreased by -0.13 ± 0.06 cm per year, or 3.64 cm over the 28 years of study. During the active period, relative humidity increased by 0.08 ± 0.03 % points per year (LM, t = 2.818, p = 0.006). Moreover, the first ordinal day of complete snowmelt increased through time (delay of 13 days within 28 years, i.e. 0.45 ± 0.18 day/year, LM, t = 2.509, p = 0.014). The only climatic variable that appeared to remain relatively invariant over the study period were wind speed and the number of snow days during the active period (LMs, wind: t = 0.070, p = 0.944, nbdays: t = 0.207, p = 0.836). In summary, the climate tended to be warmer, more humid and with a smaller snow cover during winter. Yet, this snow cover disappeared later in the spring.
630 Climate – emergence date relationships

631 *Mean relationships between climate and emergence date*

632 Among the different climate variables considered during winter (inactive phase) or the
633 previous summer (active phase), only the air temperature during previous summer had a VIF
634 greater than three (between 5.244 and 5.460 according to the sex-age category considered)
635 and was removed from the model (probably due to its high correlation with soil temperature;
636 see Appendix S1, Figure S7). The ordinal day of snowmelt and wind during summer were
637 positively related with annual mean emergence date (Fig. 5). Emergence dates were delayed
638 with a delay in the first day of snow melt for all individual categories, apart from 2-year-old
639 individuals (LMMs [13], yearlings: t = 2.211, p = 0.036, n = 84 observations, 2-year-old
640 females: t = 1.713, p = 0.099, n = 82 observations, 2-year-old males: t = 0.424, p = 0.676, n =
641 77 observations, ≥3-year-old females: t = 2.994, p = 0.007, n 83 = observations, ≥3-year-old
642 males: t = 3.958, p < 0.001, n = 83 observations). Ground squirrels therefore emerged earlier
643 when (or/and where) snow melted earlier; emergence date being delayed by approximately
644 0.25 day for a day of delayed snowmelt. Similarly, the summers with higher wind speed
645 correlated with later emergences, but only for older females (≥3-year-old, LMM, t = 2.230, p
646 = 0.037, n = 83 observations). Indeed, older females emerged around 13 days later per 1 m/s
647 increase in wind speed.

648 *Year-centered approach: are spatial differences in ED related to spatial variation in climate?*

649 Once centered per year, no strong correlations were highlighted between climate
650 variables, and we included all four of them in the model (VIFs < 1.308, see Appendix S1,
651 Figure S8). Day of snowmelt was related positively to emergence date, but only in yearlings
652 and older (≥3-year-old) females; individuals emerged later when they emerged from locations
653 where snow melted later (Fig. 6; ED around 0.2 day later for a day of delayed snowmelt, ED
654 and climate variables centered by year, LMMs [14], yearlings: t = 2.253, p = 0.024, n = 1450
655 observations, N = 1450 individuals, ≥3-year-old females: t = 2.389, p = 0.017, n = 1392
656 observations, N = 559 individuals). Males of 2-year-old and ≥3-year-old seemed to emerge
657 around 2 days later in locations where the snow depth was 1cm higher during hibernation
658 (LMMs, 2-year-old males: t = 2.303, p = 0.023, n = 280 observations, N = 280 individuals,
659 ≥3-year-old males: t = 2.085, p = 0.039, n = 693 observations, N = 300 individuals, effect
660 only marginally significant after transforming the data, Box-Cox transformation, p = 0.058).
661 Soil temperature during winter tended to have an effect on the emergence date of yearlings
662 and older males, with individuals emerging earlier from warmer hibernacula (Fig. 6, LMMs,
663 yearlings: t = -1.625, p = 0.105, n = 1450 observations, N = 1450 individuals, ≥3-year-old
664 males: t = -1.498, p = 0.137, n = 693 observations, N = 300 individuals).

665 *Spatial-centered approach: are interannual differences in ED related to interannual variation*
666 *in climate?*

667 When conducting the analysis, air temperature during the previous summer had a VIF
668 greater than three (between 4.154 and 5.249) and was removed from the model (probably due
669 to its high positive correlation with soil temperature, see Appendix S1, Figure S9). Regarding
670 temporal effects of climate on emergence date, the ordinal day of snow melt was positively
671 correlated to emergence date for all individuals excepting 2-year-old males (emergence dates
672 and climate variables centered by cluster, LMMs [15], yearlings: t = 3.015, p = 0.003, n =
673 1450 observations, N = 1450 individuals; 2-year-old females: t = 2.821, p = 0.006, n = 525
674 observations, N = 525 individuals; 2-year-old males: t = 1.275, p = 0.211, n = 280
675 observations, N = 280 individuals; ≥3-year-old females: t = 3.246, p = 0.002, n = 1392
676 observations, N = 559 individuals; ≥3-year-old males: t = 3.610, p < 0.001, n = 693
677 observations, N = 300 individuals). Thus, individuals emerged later in years when snow
678 melted later (around 0.1 day later for a day of delayed of snowmelt, see Fig. 7). Similarly,
679 active periods with high wind speed correlated positively with delayed emergence dates of
680 older individuals the next year (between 13 and 17 days delay per m/s increase, LMMs, ≥3-year-old females: t = 3.357, p = 0.001, n = 1392 observations, N = 559 individuals; ≥3-year-old males: t = 2.803, p = 0.006, n = 693 observations, N = 300 individuals). Earlier emergence date of 2-year-old males correlated with higher summer humidity (LMMs, 2-year-old males: t = -2.197, p = 0.036, n = 280 observations, N = 280 individuals), although this effect was only marginally significant after transforming the data (Box-Cox transformation, p = 0.054). Depth of snow during winter may have a positive effect on the ED of older individuals, and the emergence date of 2-year-old males seems to be negatively associated with the hibernation soil temperature (ED around 3 days earlier per degree increase of soil temperature at one-meter depth).
681
682
683
684
685
686
687
688
689

690 DISCUSSION

691 Consistent with previous studies (Dobson et al., 1992; Murie & Harris, 1982), ground
692 squirrel age and sex accounted for a substantial amount of variation (29.4%) in ED. Older
693 males emerged first from hibernation, followed by older females, 2-year-old individuals, and
694 finally yearlings. Individual identity had a significant effect in the overall pattern of
695 emergence date. This is consistent with previous knowledge on the heritability of emergence
696 date of Columbian ground squirrels (h² = 0.22 in females and 0.34 in males, Lane et al. 2011).
697 The magnitude of the individual effect (14.6% of variance) was similar to the temporal
698 variance in ED (16.1% of variance), in spite of the strong individual plasticity in ED from
699 year to year (Lane et al. 2012). These results confirm the importance of including among-
700 individuals heterogeneity in understanding population dynamics (Hamel et al. 2018).
701
702
703

ED varied over space and time (21.7% of variance), suggesting that this trait is plastic
and capable of responding to climate variation (see also Dobson, 1988; Dobson & Kjelgaard,
1985; Lane et al., 2019). Also, separating the temporal and spatial effects of climate on the
emergence, either by centering by cluster or by year, enabled us to compare different
mechanisms. Indeed, philopatric individuals (such as adult Columbian ground squirrels) are
adapted to the environmental conditions in which they live. Oppositely, interannual variations
in climate are more stochastics and might trigger plastic responses of animals. Whereas the
years effects (16.1% of variance) are consistent with previous observations that ED was
delayed through time in some populations (Lane et al., 2012; Neuhaus et al., 1999), our
results importantly show that part of the variance (5.6%) in ED was also explained by spatial
variation both within and among meadows. Such spatial heterogeneities were unlikely
explained by variations in altitude (Dobson et al. 1992), since the four study populations are
less than 70-m apart in terms of elevation. Interestingly, the variation in ED among
individuals from different locations was even higher within than among meadows, ground
squirrels from some parts of the meadows emerging earlier than others. This suggests that
local heterogeneities occur at a micro-scale (~40 m; see Results) rather than meso-scale (1-10
km). These results add to the growing literature showing that several animal and plant species
are more sensitive to microscale than wider environmental variations (Lampei et al. 2019,
Lundblad and Conway 2021, Wang and Hou 2021). Concurrently, spatial microclimate
differences were found between clusters, suggesting a potential link between ED and local
habitat conditions.

722 **Emergence date relationship to climate**

Our analyses of microclimate effects on ED focused on two distinct temporal phases:
the climate conditions during the previous summer (active period) and those during
hibernation (inactive period). First, our results showed the importance of carry-over effects of
the period preceding hibernation on ground squirrel emergence dates the subsequent year. In
particular, we found that delayed emergence occurred in older individuals when the preceding
summer was windier and associated with lower humidity (see Appendix S1, Figure S9 for
temporal correlations between climate variables). Carry-over climatic effects on animal phenology have mainly been documented in birds (Finch et al., 2014; but see Ockendon et al., 2013), fish (Wilson et al. 2021), and insects (Dingemanse and Kalkman 2008, McCauley et al. 2018). Yet, carry-over effects are widespread, including in mammals, and are thought to be primarily mediated through macronutrient supply (reviewed by Harrison et al. 2011). Climate directly affects food availability/quality or/and foraging patterns (Levy et al., 2016), which may have direct consequences on individual mass gain before hibernation. In particular, strong wind conditions may have decreased the time spent foraging for ground squirrels (Fairbanks and Dobson 2007), while low humidity associated with decreased precipitation may have resulted in reduced primary productivity and inadequate forage quantity/quality, also affecting the accumulation of fat stores (Dobson and Kjelgaard 1985, Neuhaus et al. 1999, Dobson et al. 2016). In turn, lower fat stores might have resulted in longer hibernation to avoid emerging when resources were scarce, explaining tardier emergence dates. Whereas our study suggests climatic carry-over effects on phenology, such effects are largely understudied, and yet likely widespread. For instance carry-over effects of food availability on the relationship between temperature and hibernation emergence were shown in edible dormice (Fietz et al., 2020).

Similarly, climate conditions during hibernation had marked effects on the timing of emergence of ground squirrels, either by affecting the quality of hibernation and animal energy expenditure, or acting as cues indicating the end of the winter period. Years with later snow melt resulted in later emergences for one-year-old individuals, 2-year-old females, and ≥3-year-old individuals. Similarly, spatial clusters where snow melted later resulted in yearlings and ≥3-year-old females emerging later than in other areas. This was true both at the spatial and temporal scales, and for almost all individual categories, suggesting the necessity for snow to melt and uncover growing vegetation resources before ground squirrels
754 emerge. In addition, a melting and wet snow pack makes it hard for the animals to easily
755 move around the meadow (VAV, FSD, PN; personal observations) and likely increases the
756 costs of thermoregulation, increasing overall energy expenditure at a time where their
757 remaining fat stores are likely depleted. Local climate may thus act as a strong cue on the
758 timing of ground squirrel emergence, though the question of the immediate mechanism that
759 stimulates emergence remains for future research.

760 We found a positive spatial correlation between snow depth and emergence of 2-year-
761 old and ≥3-year-old males, suggesting that males of 2-year-old and older emerged later in
762 locations where snow depth was higher during the inactive period. Higher snowpack or
763 delayed snowmelt during the inactive period may also insulate and stabilize soil temperatures
764 and hibernation conditions, decreasing fat reserve consumption, and thus allowing a longer
765 hibernation or emergence in better condition (the latter important for reproduction; Rubach et
766 al., 2016). Snowpack thickness is an important variable influencing hibernating species, such
767 as alpine marmots (Canale et al. 2016), hoary marmots (Patil et al. 2013), or northern Idaho
768 ground squirrels (Goldberg and Conway 2021). A negative correlation between the soil
769 temperature during inactive period (at one-meter depth) and the emergence date of yearlings
770 and ≥3-year-old males also seemed to apply. As for northern Idaho ground squirrels
771 (Goldberg and Conway 2021), this result possibly indicates that unfavorable hibernating
772 conditions lead to earlier emergence. For example, elevated soil temperatures prevent animals
773 from lowering their body temperature sufficiently for efficient energy conservation during
774 hibernation. In short, fat-storing and hibernating species might be sensitive to local
775 environmental conditions and especially elevated temperatures, which are directly influencing
776 energy expenditure during this period of fast.

777 Climatic variables correlating with ED differed between sex-age categories.
778 Explanations for such diverging responses may be both ecological and physiological. Young
779 (1990a) reported the existence of long post-hibernation euthermic phases only in adult males,
780 remaining several days underground before emerging despite favorable aboveground
781 conditions, possibly to allow for gonadal maturation (Williams et al. 2017). Physiological and
782 phenological differences in sex and age categories might therefore reflect distinct constraints
783 between reproductively mature individuals and non-mature individuals. Yearlings and 2-year-
784 olds that emerge later in the spring likely do so when resources become more abundant and
785 perhaps to avoid predator exposure or the social aggression related to reproduction (Murie
786 and Harris 1982, Constant et al. 2020). By contrast, mature individuals may optimize their
787 emergence date according to the trade-off between energy resources (remaining fat stores,
788 food availability) and the advantages of an early reproduction (Neuhaus 2000). Furthermore,
789 reproductive strategies might differ between males and females, with males taking advantage
790 of earlier hibernation emergence to seek more breeding opportunities (Richardson's ground
791 squirrels, Michener, 1983) and establish their territories (Manno and Dobson 2008). Thus,
792 environmental influence on a population might differ importantly between individuals with
793 different life history constraints (e.g. age classes; Bonamour et al., 2020; Pardo et al., 2013).

794 **Long-term changes in climate in relationship to emergence date**

795 Previous studies have reported shifts in animal phenology linked to changes in
796 climate, species especially responding to temperature (more so at high latitudes) and
797 precipitations (Cohen et al. 2018). In the context of climate change, advancements in spring
798 phenology have been found in the vast majority of cases (Parmesan & Yohe, 2003; but see
799 Radchuk et al., 2019). However, in studying long-term effects of climate change on animal
800 phenology, mammals are poorly represented (Parmesan & Yohe, 2003; but see Radchuck et
801 al. 2019). Among those studied, *Marmotini* species appear to respond differently to
802 directional changes in climate. The phenology of Uinta ground squirrels did not change over
803 time, despite a locally changing climate (Falvo et al. 2019). Yellow-bellied marmot have been
found to emerge earlier in response to warming air temperatures (Inouye et al. 2000), while
Columbian ground squirrels had been shown to display a trend toward later emergence,
related to late spring snowstorms (Lane et al., 2012).

Interestingly, our results are only partially consistent with previous findings on
Columbian ground squirrels. We found evidence for delayed snowmelt (13 days) over the 28
years of study despite a general warming trend (air temperature increased by 1.40°C during
the active season), and a positive effect of the first day without snow on emergence date
(except for 2-year-old individuals). Despite such a climatic trend, different phenological
responses were observed in Columbian ground squirrels of different age categories and
meadows. As an example, on Meadows C and B, all individuals (except older males on MB)
delayed emergences by 4 to 11 days over the 18 to 28 year period. Yet, no changes in
emergences were observed on DOT for any age category, and long-term temporal trends in
ED differed substantially depending on individual age and sex on Meadow A. Two-year-old
individuals advanced their emergence by 2 days over 22 years, whereas older individuals
delayed their emergence by 3 days in the same period, and no change in ED was found for
yearlings. Differences in ED trends between closely located meadows may be related to the
spatial heterogeneity of field, slope, or sun and wind exposure, with different effects of
microclimates. Similar results were found in a recent study on the North American deer
mouse (*Peromyscus maniculatus*), where temperature and photoperiod were positive cues of
breeding phenology, the relationship varying among regions (McLean and Guralnick 2021).
Whether variation in ground squirrel phenology is linked to behavioral plasticity versus long-
term adaptation remains to be seen. The small spatial scales tested here, and the diverging
responses among nearby populations might suggest phenological plasticity, as recently
suggested (Lane et al. 2019).
Many species responses to climate are likely to occur at micro-scale, rather than meso-scale. As macro-scale (e.g. sockeye salmon; *Onchorhynchus nerka*; Martins et al., 2012). As highlighted here, the relationships between phenology and climate variables are likely to differ between neighboring populations of the same species, in response to contrasting microclimate conditions, even over small spatial scales. Thus, studies focusing on climate-phenology relationships in animals would benefit by integrating fine local scale information on habitat.

## Conclusion and further perspectives

Taken together, our results point to important differences in microclimate effects on animal phenology in neighboring populations. Of course, phenology may be linked to individual fitness, and thereby the dynamics of the populations. Thus, we might also expect variation in these variables over small spatial scales. Integrating microclimatic variation in the study of climate change is likely to provide more accurate predictions of a variety of animal responses to a changing world.

## ACKNOWLEDGEMENTS

We are grateful to Edward A. Johnson (Director of the Biogeosciences Institute, University of Calgary), Judy Mappin-Buchannan and Adrienne Cunnings (Managers, Kananaskis Field Stations), and Kathreen Ruckstuhl (faculty member responsible for the R.B. Miller Field Station) for housing and facilities during fieldwork. The fieldwork was aided by many volunteers and students over the years, and we thank them for their excellent efforts.

We are specifically indebted to Jan O. Murie for initiating the long-term research on Columbian ground squirrels, and to David A. Boag for his earlier work on the species. Jan O. Murie provided insightful comments on the manuscript, and we are grateful for his continued
advice, help and friendship over the years. This research was supported by a USA National

Science Foundation grant (DEB-0089473) and a fellowship grant from the Institute of

Advanced Studies of the University of Strasbourg to F.S.D., a Natural Sciences and

Engineering Research Council of Canada grant to J.O. Murie, a CNRS Projet International de

Coopération Scientifique grant (PICS-07143), an AXA Postdoctoral Research Fellowship, a

Fyssen Foundation Postdoctoral Fellowship, and a Fyssen Foundation research grant to

V.A.V., and an Initiative d'Excellence (IDEX attractivité) research grant to C. Saraux. F.S.D.

thanks the Région Grand Est and the Eurométropole de Strasbourg for the award of a

Gutenberg Excellence Chair during the time of writing. A.T. was supported by a PhD

scholarship from the Ministère de l'Enseignement Supérieur, de la Recherche et de

l'Innovation.

## 863 ETHICAL APPROVAL

All applicable international, national, and/or institutional guidelines for the care and

use of animals were followed. All procedures carried out in the field were approved by

Auburn University and the University of Calgary (IACUC protocols) and the University of

Calgary. Authorization for conducting research and collecting samples in the Sheep River

Provincial Park were obtained from Alberta Environment and Parks and Alberta Fish &

Wildlife.

## 870 REFERENCES

Alberta Parks. 2008. Sheep River Provincial Park and Bluerock Wildland Provincial Park

management plan - Open Government. https://open.alberta.ca/dataset/9780778574156.

Arnaud, C. M., F. S. Dobson, and J. O. Murie. 2012. Philopatry and within-colony

movements in Columbian ground squirrels. Molecular Ecology 21:493–504.

Barry, R. G. 1992. Mountain Weather and Climate. Psychology Press.
876 Bonamour, S., L.-M. Chevin, D. Réale, C. Teplitsky, and A. Charmantier. 2020. Age-

877 dependent phenological plasticity in a wild bird. Journal of Animal Ecology 89:2733–

878 2741.

879 Box, G. E. P., and D. R. Cox. 1964. An Analysis of Transformations. Journal of the Royal

880 Statistical Society: Series B (Methodological) 26:211–243.

881 Burnham, K. P., and D. R. Anderson. 2004. Multimodel Inference: Understanding AIC and

882 BIC in Model Selection. Sociological Methods & Research 33:261–304.

883 Canale, C. I., A. Ozgul, D. Allainé, and A. Cohas. 2016. Differential plasticity of size and

884 mass to environmental change in a hibernating mammal. Global Change Biology

885 22:3286–3303.

886 Charmantier, A., R. H. McCleery, L. R. Cole, C. Perrins, L. E. B. Kruuk, and B. C. Sheldon.

887 2008. Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird

888 Population. Science 320:800–803.

889 Chmura, H. E., H. M. Kharouba, J. Ashander, S. M. Ehlman, E. B. Rivest, and L. H. Yang.

890 2019. The mechanisms of phenology: the patterns and processes of phenological

891 shifts. Ecological Monographs 89:e01337.

892 Cleland, E. E., I. Chuine, A. Menzel, H. A. Mooney, and M. D. Schwartz. 2007. Shifting

893 plant phenology in response to global change. Trends in Ecology & Evolution 22:357–

894 365.

895 Cohen, J. M., M. J. Lajeunesse, and J. R. Rohr. 2018. A global synthesis of animal

896 phenological responses to climate change. Nature Climate Change 8:224–228.

897 Constant, T., S. Giroud, V. A. Viblanc, M. L. Tissier, P. Bergeron, F. S. Dobson, and C.

898 Habold. 2020. Integrating Mortality Risk and the Adaptiveness of Hibernation.

899 Frontiers in Physiology 11:706.
900 Davis, D. E. 1976. Hibernation and Circannual Rhythms of Food Consumption in Marmots
901 and Ground Squirrels. The Quarterly Review of Biology 51:477–514.

902 Delignette-Muller, M. L., and C. Dutang. 2015. fitdistrplus: An R package for fitting
903 distributions. Journal of statistical software 64:1–34.

904 Dingemanse, N. J., and V. J. Kalkman. 2008. Changing temperature regimes have advanced
905 the phenology of Odonata in the Netherlands. Ecological Entomology 33:394–402.

906 Dobson, F. S. 1988. The limits of phenotypic plasticity in life histories of Columbian ground
907 squirrels. Pages 193–210 Evolution of life histories of mammals.

908 Dobson, F. S., M. J. Badry, and C. Geddes. 1992. Seasonal activity and body mass of
909 Columbian ground squirrels. Canadian Journal of Zoology 70:1364–1368.

910 Dobson, F. S., and J. D. Kjelgaard. 1985. The influence of food resources on life history in
911 Columbian ground squirrels. Canadian Journal of Zoology 63:2105–2109.

912 Dobson, F. S., J. E. Lane, M. Low, and J. O. Murie. 2016. Fitness implications of seasonal
913 climate variation in Columbian ground squirrels. Ecology and Evolution 6:5614–5622.

914 Durant, J. M., D. Ø. Hjermann, G. Ottersen, and N. C. Stenseth. 2007. Climate and the match
915 or mismatch between predator requirements and resource availability. Climate
916 Research 33:271–283.

917 Fairbanks, B., and F. S. Dobson. 2007. Mechanisms of the group-size effect on vigilance in
918 Columbian ground squirrels: dilution versus detection. Animal Behaviour 73:115–123.

919 Falvo, C. A., D. N. Koons, and L. M. Aubry. 2019. Seasonal climate effects on the survival of
920 a hibernating mammal. Ecology and Evolution 9:3756–3769.

921 Festa-Bianchet, M., and D. A. Boag. 1982. Territoriality in adult female Columbian ground
922 squirrels. Canadian Journal of Zoology 60:1060–1066.

923 Fietz, J., F. Langer, and W. Schlund. 2020. They like it cold, but only in winter: Climate-
924 mediated effects on a hibernator. Functional Ecology 34:2098–2109.
Finch, T., J. W. Pearce-Higgins, D. I. Leech, and K. L. Evans. 2014. Carry-over effects from passage regions are more important than breeding climate in determining the breeding phenology and performance of three avian migrants of conservation concern. *Biodiversity and Conservation* 23:2427–2444.

Geiser, F. 2013. Hibernation. Current Biology 23:R188–R193.

Goldberg, A. R., and C. J. Conway. 2021. Hibernation behavior of a federally threatened ground squirrel: climate change and habitat selection implications. *Journal of Mammalogy* 102:574–587.

Graae, B. J., P. D. Frenne, A. Kolb, J. Brunet, O. Chabrerie, K. Verheyen, N. Pepin, T. Heinken, M. Zobel, A. Shevtsova, I. Nijs, and A. Milbau. 2012. On the use of weather data in ecological studies along altitudinal and latitudinal gradients. *Oikos* 121:3–19.

Gultepe, I. 2015. Mountain Weather: Observation and Modeling. Pages 229–312 *Advances in Geophysics*. Elsevier.

Hamel, S., J.-M. Gaillard, N. G. Yoccoz, R. D. Bassar, S. Bouwhuis, H. Caswell, M. Douhard, E. J. Gangloff, O. Gimenez, P. C. Lee, I. M. Smallegange, U. K. Steiner, O. Vedder, and Y. Vindenes. 2018. General conclusion to the special issue *Moving forward on individual heterogeneity*. *Oikos* 127:750–756.

Harris, M. A., and J. O. Murie. 1984. Inheritance of Nest Sites in Female Columbian Ground Squirrels on JSTOR. *Behavioral Ecology and Sociobiology* 15:97–102.

Harrison, X. A., J. D. Blount, R. Inger, D. R. Norris, and S. Bearhop. 2011. Carry-over effects as drivers of fitness differences in animals. *Journal of Animal Ecology* 80:4–18.

Hengl, T., J. M. de Jesus, G. B. M. Heuvelink, M. R. Gonzalez, M. Kilibarda, A. Blagotić, W. Shangguan, M. N. Wright, X. Geng, B. Bauer-Marschallinger, M. A. Guevara, R. Vargas, R. A. MacMillan, N. H. Batjes, J. G. B. Leenaars, E. Ribeiro, I. Wheeler, S.
949 Mantel, and B. Kempen. 2017. SoilGrids250m: Global gridded soil information based
950 on machine learning. PLOS ONE 12:e0169748.

951 Hollister J, Shah T, Robitaille A, Beck M, and Johnson M. 2020. elevatr: Access Elevation
952 Data from Various APIs. R package version 0.3.1.

953 Huggett, R. J. 1995. Geoecology: An Evolutionary Approach. Taylor & Francis US.

954 Humphries, M. M., D. W. Thomas, and D. L. Kramer. 2003. The Role of Energy Availability
955 in Mammalian Hibernation: A Cost-Benefit Approach. Physiological and Biochemical
956 Zoology 76:165-179.

957 Humphries, M. M., D. W. Thomas, and J. R. Speakman. 2002. Climate-mediated energetic
958 constraints on the distribution of hibernating mammals. Nature 418:313-316.

959 Humphries, M. M., J. Umbanhowar, and K. S. McCann. 2004. Bioenergetic Prediction of
960 Climate Change Impacts on Northern Mammals1. Integrative and Comparative
961 Biology 44:152-162.

962 Inouye, D. W., B. Barr, K. B. Armitage, and B. D. Inouye. 2000. Climate change is affecting
963 altitudinal migrants and hibernating species. Proceedings of the National Academy of
964 Sciences 97:1630-1633.

965 James, F. C., R. F. Johnston, N. O. Wamer, G. J. Niemi, and W. J. Boecklen. 1984. The
966 Grinnellian Niche of the Wood Thrush. The American Naturalist 124:17-47.

967 Kalnay, E., M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell, S. Saha,
968 G. White, J. Woollen, Y. Zhu, M. Chelliah, W. Ebisuzaki, W. Higgins, J. Janowiak, K.
969 C. Mo, C. Ropelewski, J. Wang, A. Leetmaa, R. Reynolds, R. Jenne, and D. Joseph.
970 1996. The NCEP/NCAR 40-Year Reanalysis Project. Bulletin of the American
971 Meteorological Society 77:437-472.

972 Kearney, M. R. 2020. How will snow alter exposure of organisms to cold stress under climate
973 warming? Global Ecology and Biogeography 29:1246-1256.
974 Kearney, M. R., P. K. Gillingham, I. Bramer, J. P. Duffy, and I. M. D. Maclean. 2020. A
975 method for computing hourly, historical, terrain-corrected microclimate anywhere on
976 earth. Methods in Ecology and Evolution 11:38–43.

977 Kreyling, J., and H. A. L. Henry. 2011. Vanishing winters in Germany: soil frost dynamics
978 and snow cover trends, and ecological implications. Climate Research 46:269–276.

979 Kunkel, K. E., D. R. Easterling, K. Hubbard, and K. Redmond. 2004. Temporal variations in
980 frost-free season in the United States: 1895–2000. Geophysical Research Letters 31.

981 Lampei, C., J. Wunder, T. Wilhalm, and K. J. Schmid. 2019. Microclimate predicts frost
982 hardiness of alpine Arabidopsis thaliana populations better than elevation. Ecology
983 and Evolution 9:13017–13029.

984 Lane, J. E., Z. J. Czenze, R. Findlay-Robinson, and E. Bayne. 2019. Phenotypic Plasticity and
985 Local Adaptation in a Wild Hibernator Evaluated through Reciprocal Translocation.
986 The American Naturalist 197:516–528.

987 Lane, J. E., L. E. B. Kruuk, A. Charmantier, J. O. Murie, D. W. Coltman, M. Buoro, S.
988 Raveh, and F. S. Dobson. 2011. A quantitative genetic analysis of hibernation
989 emergence date in a wild population of Columbian ground squirrels. Journal of
990 Evolutionary Biology 24:1949–1959.

991 Lane, J. E., L. E. B. Kruuk, A. Charmantier, J. O. Murie, and F. S. Dobson. 2012. Delayed
992 phenology and reduced fitness associated with climate change in a wild hibernator.
993 Nature 489:554–557.

994 Lê, S., J. Josse, and F. Husson. 2008. FactoMineR: an R package for multivariate analysis.
995 Journal of statistical software 25:1–18.

996 Lembrechts, J. J., and J. Lenoir. 2020. Microclimatic conditions anywhere at any time! Global
997 Change Biology 26:337–339.
998

Lenth, R. V. 2016. Least-Squares Means: The R Package lsmeans. Journal of Statistical

999

Software 69:1–33.

1000

Levy, O., T. Dayan, W. P. Porter, and N. Kronfeld-Schor. 2016. Foraging Activity Pattern Is

1001

Shaped by Water Loss Rates in a Diurnal Desert Rodent. The American Naturalist

1002

188:205–218.

1003

Lundblad, C. G., and C. J. Conway. 2021. Nest microclimate and limits to egg viability

1004

explain avian life-history variation across latitudinal gradients. Ecology 102:e03338.

1005

Maclean, I. M. D. 2020. Predicting future climate at high spatial and temporal resolution.

1006

Global Change Biology 26:1003–1011.

1007

Maclean, I. M. D., A. J. Suggitt, R. J. Wilson, J. P. Duffy, and J. J. Bennie. 2017. Fine-scale

1008

climate change: modelling spatial variation in biologically meaningful rates of

1009

warming. Global Change Biology 23:256–268.

1010

Manno, T. G., and F. S. Dobson. 2008. Why are Male Columbian Ground Squirrels

1011

Territorial? Ethology 114:1049–1060.

1012

Martay, B., M. J. Brewer, D. A. Elston, J. R. Bell, R. Harrington, T. M. Brereton, K. E.

1013

Barlow, M. S. Botham, and J. W. Pearce-Higgins. 2017. Impacts of climate change on

1014

national biodiversity population trends. Ecography 40:1139–1151.

1015

Martins, E. G., S. G. Hinch, S. J. Cooke, and D. A. Patterson. 2012. Climate effects on

1016

growth, phenology, and survival of sockeye salmon (Oncorhynchus nerka): a synthesis

1017

of the current state of knowledge and future research directions. Reviews in Fish

1018

Biology and Fisheries 22:887–914.

1019

Matheron, G. 1963. Principles of geostatistics. Economic Geology 58:1246–1266.

1020

McCafferty, D. J., G. Pandraud, J. Gilles, M. Fabra-Puchol, and P.-Y. Henry. 2017. Animal

1021

thermoregulation: a review of insulation, physiology and behaviour relevant to

1022

temperature control in buildings. Bioinspiration & Biomimetics 13:011001.
1023 McCauley, S. J., J. I. Hammond, and K. E. Mabry. 2018. Simulated climate change increases larval mortality, alters phenology, and affects flight morphology of a dragonfly. Ecosphere 9:e02151.

1026 McGuire, L. P., E. M. Johnson, W. F. Frick, and J. G. Boyles. 2021. Temperature alone is insufficient to understand hibernation energetics. The Journal of Experimental Biology.

1029 McLean, B. S., and R. P. Guralnick. 2021. Digital biodiversity data sets reveal breeding phenology and its drivers in a widespread North American mammal. Ecology 102:e03258.

1032 Menzel, A. 2002. Phenology: Its Importance to the Global Change Community. Climatic Change 54:379–385.

1034 Menzel, A., and P. Fabian. 1999. Growing season extended in Europe. Nature 397:659–659.

1035 Michener, G. R. 1983. Spring emergence schedules and vernal behavior of Richardson's ground squirrels: why do males emerge from hibernation before females? Behavioral Ecology and Sociobiology 14:29–38.

1038 Moore, M. P., and R. A. Martin. 2019. On the evolution of carry-over effects. Journal of Animal Ecology 88:1832–1844.

1040 Murie, J. O., and D. A. Boag. 1984. The Relationship of Body Weight to Overwinter Survival in Columbian Ground Squirrels. Journal of Mammalogy 65:688–690.

1042 Murie, J. O., and M. A. Harris. 1978. Territoriality and dominance in male Columbian ground squirrels (Spermophilus columbianus). Canadian Journal of Zoology 56:2402–2412.

1044 Murie, J. O., and M. A. Harris. 1982. Annual Variation of Spring Emergence and Breeding in Columbian Ground Squirrels (Spermophilus columbianus). Journal of Mammalogy 63:431–439.
1047 Murie, J. O., and M. A. Harris. 2011. Social interactions and dominance relationships

1048 between female and male Columbian ground squirrels. Canadian Journal of Zoology

1049 66:1414–1420.

1050 Nagy, L., and G. Grabherr. 2009. The Biology of Alpine Habitats. OUP Oxford.

1051 Neuhaus, P. 2000. Timing of Hibernation and Molt in Female Columbian Ground Squirrels.

1052 Journal of Mammalogy 81:571–577.

1053 Neuhaus, P., R. Bennett, and A. Hubbs. 1999. Effects of a late snowstorm and rain on survival

1054 and reproductive success in Columbian ground squirrels (Spermophilus columbianus).

1055 Canadian Journal of Zoology 77:879–884.

1056 Norquay, K. J. O., and C. K. R. Willis. 2014. Hibernation phenology of Myotis lucifugus.

1057 Journal of Zoology 294:85–92.

1058 Norris, D. R. 2005. Carry-over effects and habitat quality in migratory populations. Oikos

1059 109:178–186.

1060 Ockendon, N., D. Leech, and J. W. Pearce-Higgins. 2013. Climatic effects on breeding

1061 grounds are more important drivers of breeding phenology in migrant birds than carry-

1062 over effects from wintering grounds. Biology Letters 9:20130669.

1063 O'Rourke and Hatcher. 2013. A Step-by-Step Approach to Using SAS for Factor Analysis

1064 and Structural Equation Modeling, Second Edition. SAS Institute.

1065 Ozgul, A., D. Z. Childs, M. K. Oli, K. B. Armitage, D. T. Blumstein, L. E. Olson, S.

1066 Tuljapurkar, and T. Coulson. 2010. Coupled dynamics of body mass and population

1067 growth in response to environmental change. Nature 466:482–485.

1068 Pachauri, R. K., M. R. Allen, V. R. Barros, J. Broome, W. Cramer, R. Christ, J. A. Church, L.

1069 Clarke, Q. Dahe, P. Dasgupta, N. K. Dubash, O. Edenhofer, I. Elgizouli, C. B. Field,

1070 P. Forster, P. Friedlingstein, J. Fuglestvedt, L. Gomez-Echeverri, S. Hallegatte, G.

1071 Hegerl, M. Howden, K. Jiang, B. Jimenez Cisneroz, V. Kattsov, H. Lee, K. J. Mach, J.
1072 Marotzke, M. D. Mastrandrea, L. Meyer, J. Minx, Y. Mulugetta, K. O'Brien, M.

1073 Oppenheimer, J. J. Pereira, R. Pichs-Madruga, G.-K. Plattner, H.-O. Pörtner, S. B.

1074 Power, B. Preston, N. H. Ravindranath, A. Reisinger, K. Riahi, M. Rusticucci, R.

1075 Scholes, K. Seyboth, Y. Sokona, R. Stavins, T. F. Stocker, P. Tschakert, D. van Vuuren, and J.-P. van Ypserle. 2014b. Climate Change 2014: Synthesis Report.

1077 Contribution of Working Groups I, II and III to the Fifth Assessment Report of the

1078 Intergovernmental Panel on Climate Change. Page 151 (R. K. Pachauri and L. Meyer, Eds.) EPIC3Geneva, Switzerland, IPCC, 151 p., pp. 151, ISBN: 978-92-9169-143-2.

1080 IPCC, Geneva, Switzerland.

1081 Pardo, D., C. Barbraud, M. Authier, and H. Weimerskirch. 2013. Evidence for an age-dependent influence of environmental variations on a long-lived seabird's life-history traits. Ecology 94:208-220.

1084 Parmesan, C., and G. Yohe. 2003. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421:37-42.

1086 Paruchuri, S., A. T. Smith, Z. Fan, and F. S. Dobson. 2019. Microhabitat use by plateau pikas: living on the edge. Journal of Mammalogy 100:1221-1228.

1088 Patil, V. P., S. F. Morrison, T. J. Karels, and D. S. Hik. 2013. Winter weather versus group thermoregulation: what determines survival in hibernating mammals? Oecologia 173:139-149.

1091 Pires, A. P. F., D. S. Srivastava, N. A. C. Marino, A. A. M. MacDonald, M. P. Figueiredo-Barros, and V. F. Farjalla. 2018. Interactive effects of climate change and biodiversity loss on ecosystem functioning. Ecology 99:1203-1213.

1094 Potter, K. A., H. A. Woods, and S. Pincebourde. 2013. Microclimatic challenges in global change biology. Global Change Biology 19:2932-2939.
1096 Pretzlaff, I., and K. H. Dausmann. 2012. Impact of Climatic Variation on the Hibernation

1097 Physiology of Muscardinus avellanarius. Pages 85–97 in T. Ruf, C. Bieber, W.

1098 Arnold, and E. Millesi, editors. Living in a Seasonal World: Thermoregulatory and

1099 Metabolic Adaptations. Springer, Berlin, Heidelberg.

1100 Radchuk, V., T. Reed, C. Teplitsky, M. van de Pol, A. Charmantier, C. Hassall, P. Adamík, F.

1101 Adriaensen, M. P. Ahola, P. Arcese, J. Miguel Avilés, J. Balbontin, K. S. Berg, A.

1102 Borras, S. Burthe, J. Clobert, N. Dehnhard, F. de Lope, A. A. Dhondt, N. J.

1103 Dingemanse, H. Doi, T. Eeva, J. Fickel, I. Filella, F. Fossøy, A. E. Goodenough, S. J.

1104 G. Hall, B. Hansson, M. Harris, D. Hasselquist, T. Hickler, J. Joshi, H. Kharouba, J.

1105 G. Martínez, J.-B. Mihoub, J. A. Mills, M. Molina-Morales, A. Moksnes, A. Ozgul, D.

1106 Parejo, P. Pilard, M. Poisbleau, F. Rousset, M.-O. Rödel, D. Scott, J. C. Senar, C.

1107 Stefanescu, B. G. Stokke, T. Kusano, M. Tarka, C. E. Tarwater, K. Thonicke, J.

1108 Thorley, A. Wilting, P. Tryjanowski, J. Merilä, B. C. Sheldon, A. Pape Møller, E.

1109 Matthysen, F. Janzen, F. S. Dobson, M. E. Visser, S. R. Beissinger, A. Courtiol, and S.

1110 Kramer-Schadt. 2019. Adaptive responses of animals to climate change are most
1111 likely insufficient. Nature Communications 10:1–14.

1112 Rivoirard, J., J. Simmonds, K. G. Foote, P. Fernandes, and N. Bez. 2008. Geostatistics for
1113 Estimating Fish Abundance. John Wiley & Sons.

1114 Root, T. L., J. T. Price, K. R. Hall, S. H. Schneider, C. Rosenzweig, and J. A. Pounds. 2003.
1115 Fingerprints of global warming on wild animals and plants. Nature 421:57–60.

1116 Rubach, K., M. Wu, A. Abebe, F. S. Dobson, J. O. Murie, and V. A. Viblanc. 2016. Testing
1117 the reproductive and somatic trade-off in female Columbian ground squirrels. Ecology
1118 and Evolution 6:7586–7595.
1119 Saino, N., R. Ambrosini, M. Caprioli, A. Romano, M. Romano, D. Rubolini, C. Scandolara, and F. Liechti. 2017. Sex-dependent carry-over effects on timing of reproduction and fecundity of a migratory bird. Journal of Animal Ecology 86:239-249.

1122 Saleeba, K., B. L. Phillips, M. O'Shea, and M. R. Kearney. 2020. Using Biophysical Models to Improve Survey Efficiency for Cryptic Ectotherms. The Journal of Wildlife Management 84:1185-1195.

1125 Saraux, C., J.-M. Fromentin, J.-L. Bigot, J.-H. Bourdeix, M. Morfin, D. Roos, E. V. Beveren, and N. Bez. 2014. Spatial Structure and Distribution of Small Pelagic Fish in the Northwestern Mediterranean Sea. PLOS ONE 9:e111211.

1128 Shaw, W. T. 1926. The Storing Habit of the Columbian Ground Squirrel. The American Naturalist 60:367-373.

1130 Sparks, T. H., and T. J. Yates. 1997. The effect of spring temperature on the appearance dates of British butterflies 1883-1993. Ecography 20:368-374.

1132 Strode, P. K. 2003. Implications of climate change for North American wood warblers (Parulidae). Global Change Biology 9:1137-1144.

1134 Suggitt, A. J., P. K. Gillingham, J. K. Hill, B. Huntley, W. E. Kunin, D. B. Roy, and C. D. Thomas. 2011. Habitat microclimates drive fine-scale variation in extreme temperatures. Oikos 120:1-8.

1137 Varner, J., and M. D. Dearing. 2014. The Importance of Biologically Relevant Microclimates in Habitat Suitability Assessments. PLoS ONE 9:e104648.

1139 Visser, M. E., and C. Both. 2005. Shifts in phenology due to global climate change: the need for a yardstick. Proceedings of the Royal Society B: Biological Sciences 272:2561-2569.
1142 Visser, M. E., and L. J. M. Holleman. 2001. Warmer springs disrupt the synchrony of oak and winter moth phenology. Proceedings of the Royal Society of London. Series B: Biological Sciences 268:289–294.

1145 Walther, G.-R., E. Post, P. Convey, A. Menzel, C. Parmesan, T. J. C. Beebee, J.-M. Fromentin, O. Hoegh-Guldberg, and F. Bairlein. 2002. Ecological responses to recent climate change. Nature 416:389–395.

1148 Wang, S., and F. Hou. 2021. Burrow characteristics and ecological significance of Marmota himalayana in the northeastern Qinghai-Tibetan Plateau. Ecology and Evolution 11:9100–9109.

1151 Wiggett, D. R., D. A. Boag, and A. D. R. Wiggett. 1989. Movements of intercolony natal dispersers in the Columbian ground squirrel. Canadian Journal of Zoology 67:1447–1452.

1154 Williams, C. T., B. M. Barnes, G. J. Kenagy, and C. L. Buck. 2014. Phenology of hibernation and reproduction in ground squirrels: integration of environmental cues with endogenous programming. Journal of Zoology 292:112–124.

1157 Williams, C. T., C. L. Buck, M. J. Sheriff, M. M. Richter, J. S. Krause, and B. M. Barnes. 2017. Sex-Dependent Phenological Plasticity in an Arctic Hibernator. American Naturalist 190:854–859.

1160 Wilson, S. M., T. W. Buehrens, J. L. Fisher, K. L. Wilson, and J. W. Moore. 2021. Phenological mismatch, carryover effects, and marine survival in a wild steelhead trout Oncorhynchus mykiss population. Progress in Oceanography 193:102533.

1163 Winder, M., and D. E. Schindler. 2004. Climate Change Uncouples Trophic Interactions in an Aquatic Ecosystem. Ecology 85:2100–2106.

1165 WMO. 2015, November 26. Climate. https://public.wmo.int/en/our-mandate/climate.
1166 Woillez, M., J.-C. Poulard, J. Rivoirard, P. Petitgas, and N. Bez. 2007. Indices for capturing
spatial patterns and their evolution in time, with application to European hake
1167
1168 (Merluccius merluccius) in the Bay of Biscay. ICES Journal of Marine Science
1169 64:537–550.

1170 Wood, S. 2017. Generalized Additive Models: An Introduction with R. Chapman and
1171 Hall/CRC.

1172 Young, P. J. 1990a. Structure, Location and Availability of Hibernacula of Columbian
1173 Ground Squirrels (Spermophilus columbianus). The American Midland Naturalist
1174 123:357–364.

1175 Young, P. J. 1990b. Hibernating patterns of free-ranging Columbian ground squirrels.
1176 Oecologia 83:504–511.

1177 Zellweger, F., P. De Frenne, J. Lenoir, D. Rocchini, and D. Coomes. 2019. Advances in
1178 Microclimate Ecology Arising from Remote Sensing. Trends in Ecology & Evolution
1179 34:327–341.

1180 Zuur, A. F., E. N. Ieno, and C. S. Elphick. 2010. A protocol for data exploration to avoid
1181 common statistical problems. Methods in Ecology and Evolution 1:3–14.

1182
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Study site composed of meadows and clusters within meadows. (A) The four different meadows where populations were monitored (grey areas on the map). (B) Emergence locations (black dots) grouped into spatial clusters on each meadow. Clusters are represented by minimum convex polygons and colored according to their mean emergence date. Lighter colors correspond to earlier mean emergence date of the cluster (yellow and light orange) and darker colors correspond to later mean emergence date of the cluster (dark orange and red). Mean ordinal emergence date on clusters ranged from 118 to 127 on MA, from 118 to 130 on MB, from 120 to 124 on MC and from 114 to 124 on DOT.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Interannual variability in emergence date. Temporal patterns of emergence dates (in ordinal day) on each meadow (MA, MB, MC, DOT) and for different individual categories. Yearlings (green squares, 1), 2-year-old individuals (orange circles, 2), older females (blue triangles, 3+ Female), and older males (pink diamonds, 3+ Male) are represented. Plotted lines correspond to significant temporal trends and their 95% confidence interval according to GAMs [5], excepting the temporal pattern of older males on MB and MC which is non-significant (*ns*).</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3. Spatial variation of climate conditions.** Left side: Climate variables, considered either during the active or inactive period, included in the Principal Component Analysis along dimensions 1 and 2 (A) and dimensions 1 and 3 (C). The number of days with snowfalls during the active period (Days_snow) and the ordinal day of snowmelt (Day_snowmelt) were not used to estimate the PCA but projected on it. Right side: Annual mean values of climate conditions for each cluster of the four meadows. Clusters are represented by minimum convex polygons colored by meadows (MA in shades of blue, MB in shades of green, MC in shades of red and DOT in shades of yellow), with their center of gravity represented by larger dots.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4. Climate interannual variability.** Annual mean values of climate variables averaged by meadow. The lines correspond to temporal patterns of climate variables estimated by the GAM [9] on the four meadows (MA in blue, MB in green, MC in purple, and DOT in red) between 1992 and 2019. Climate time series were divided into two periods: inactive from the 1st of August to the 15th of April (right column) and then active the rest of the year (left column) and averaged by year. Where no difference was found between meadows (all variables except wind and humidity), a single prediction line was plotted for all locations (in black) with its 95% confidence interval (shaded area).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 5. Climate effects on hibernation emergence date of different individual categories.</figcaption></figure>

From (A) to (E): Estimates and 95% confidence intervals from the full LMM [13] explaining annual mean emergence dates per meadow are presented. Significant effects are presented in orange while non-significant ones are blue. (F) Annual mean emergence date as a function of the wind speed during the previous active period for ≥3-year-old females and model predictions (predicted line and 95% CI around). From (G) to (I): Annual mean emergence date as a function of the ordinal day of snow melt and model predictions for one-year-old individuals (G), ≥3-year-old females (H) and ≥3-year-old males (I).
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 6.** Spatial climate effects on centered emergence dates of different individual categories. From (A) to (E): Estimates and 95% confidence intervals from the full LMMs [14] are presented with significant effects in orange. From (F) to (I): Spatial effect of climate on emergence date (year-centered observations to remove temporal variance). Model predictions are indicated along with their 95% confidence intervals (shaded area) on top of observations. Significant spatial effect of the ordinal day of snowmelt on emergence date of (F) yearlings and (H) older females (≥3-year-old and more). Significant spatial effect of the winter snow depth on the emergence date of (G) 2-year-old males and (I) ≥3-year-old males.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 7.** Temporal climate effects on centered emergence dates of different individual categories. From (A) to (E): Estimates and 95% confidence intervals from the full LMMs [15] are presented with significant effects in orange. From (F) to (L): Temporal effect of climate on emergence date (all cluster-centered to remove the spatial variance). Model predictions are indicated along with their 95% confidence intervals (shaded area) on top of observations. Significant temporal effects of the ordinal day of snowmelt on the ED of (F) yearlings, (G) 2-year-old females, (J) ≥3-year-old females, and (L) ≥3-year-old males. Significant temporal effects of (H) relative humidity during preceding summer on the ED of 2-year-old males, and of wind speed on the ED of ≥3-year-old (I) females and (K) males.</figcaption></figure>
# Appendix S1: Supplementary Figures and Tables

**Journal Name:** Ecosphere

**Article Title:** Integrating microclimatic variation in phenological responses to climate change: a 28-year study in a hibernating mammal

**Authors:** Anouch TAMIAN<fnref n="1" />ª, Vincent A VIBLANC<fnref n="1" />*, F Stephen DOBSON<fnref n="1" />,<fnref n="2" />, Peter NEUHAUS<fnref n="3" />, Tracey L HAMMER<fnref n="1" />,<fnref n="3" />, Anna P NESTEROVA<fnref n="4" />, Shirley RAVEH<fnref n="5" />, Amy L SKIBIEL<fnref n="6" />, David BROUSSARD<fnref n="7" />, Theodore G MANNO<fnref n="8" />, Nandini RAJAMANI<fnref n="9" />, Claire SARAUX<fnref n="1" />*

<footnote n="1">¹Institut Pluridisciplinaire Hubert Curien, CNRS, Département Ecologie, Physiologie et Ethologie, 23 Rue du Loess, 67037 Strasbourg, France</footnote>

<footnote n="2">²Department of Biological Sciences, Auburn University, Auburn, AL 36849, USA</footnote>

<footnote n="3">³Department of Biological Sciences, University of Calgary, Canada</footnote>

<footnote n="4">⁴FaunaStats, 16 avenue de l'Europe, Immeuble SXB1, 67300 Schiltigheim, France</footnote>

<footnote n="5">⁵Institute of Biodiversity, Animal Health and Comparative Medicine, University of Glasgow, Glasgow, UK</footnote>

<footnote n="6">⁶Department of Animal, Veterinary and Food Sciences, University of Idaho, Moscow, USA</footnote>

<footnote n="7">⁷Department of Biology, Lycoming College, Williamsport, PA 17701, USA</footnote>

<footnote n="8">⁸Science Department, Catalina Foothills High School, 4300 E. Sunrise Drive, Tucson, AZ 85718 USA</footnote>

<footnote n="9">⁹Indian Institute of Science Education and Research, Tirupati Andhra Pradesh, India</footnote>

<footnote>* Co-senior authors</footnote>

<footnote>ªCorresponding author: anouch.tamian@iphc.cnrs.fr</footnote>
<table><caption>**Table S1.** Table of emergence date observations per year on the four different meadows. The number of observations is presented for each individual (sex-age) category: 1-year-old individuals (1), 2-year-old individuals (2), ≥3-year-old females (3+Fem), and ≥3-year-old males (3+Mal).</caption><thead><tr><th rowspan="2">Year</th><th colspan="4">Meadow A</th><th colspan="4">Meadow B</th><th colspan="4">Meadow C</th><th colspan="4">Meadow DOT</th></tr><tr><th>1</th><th>2</th><th>3+Fem</th><th>3+Mal</th><th>1</th><th>2</th><th>3+Fem</th><th>3+Mal</th><th>1</th><th>2</th><th>3+Fem</th><th>3+Mal</th><th>1</th><th>2</th><th>3+Fem</th><th>3+Mal</th></tr></thead><tbody><tr><td>1992</td><td></td><td></td><td></td><td></td><td>17</td><td>1</td><td>4</td><td>6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1993</td><td></td><td></td><td></td><td></td><td>12</td><td>13</td><td>4</td><td>2</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1994</td><td>74</td><td>22</td><td>65</td><td>31</td><td>22</td><td>9</td><td>11</td><td>6</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1995</td><td>73</td><td>39</td><td>66</td><td>33</td><td>33</td><td>10</td><td>15</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1996</td><td>10</td><td>43</td><td>68</td><td>32</td><td>22</td><td>13</td><td>18</td><td>9</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1997</td><td>8</td><td>9</td><td>71</td><td>22</td><td>5</td><td>5</td><td>24</td><td>9</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1998</td><td>25</td><td>6</td><td>58</td><td>24</td><td>32</td><td>3</td><td>25</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1999</td><td>23</td><td>13</td><td>43</td><td>20</td><td>43</td><td>17</td><td>26</td><td>8</td><td>17</td><td></td><td></td><td></td><td>2</td><td></td><td></td><td></td></tr><tr><td>2000</td><td>5</td><td>14</td><td>27</td><td>17</td><td>28</td><td>22</td><td>25</td><td>10</td><td>19</td><td>10</td><td></td><td></td><td>1</td><td></td><td></td><td></td></tr><tr><td>2001</td><td>3</td><td>11</td><td>26</td><td>16</td><td>35</td><td>21</td><td>39</td><td>13</td><td>26</td><td>2</td><td>8</td><td>11</td><td>5</td><td></td><td>10</td><td></td></tr><tr><td>2002</td><td></td><td>1</td><td>24</td><td>12</td><td>53</td><td>15</td><td>41</td><td>11</td><td></td><td></td><td></td><td></td><td>8</td><td>6</td><td>22</td><td>9</td></tr><tr><td>2003</td><td></td><td></td><td></td><td></td><td>4</td><td>6</td><td>21</td><td>6</td><td></td><td></td><td></td><td></td><td>18</td><td>9</td><td>32</td><td>13</td></tr><tr><td>2004</td><td></td><td></td><td></td><td></td><td>7</td><td>4</td><td>17</td><td>4</td><td>5</td><td>5</td><td>10</td><td>5</td><td>16</td><td>7</td><td>29</td><td></td></tr><tr><td>2005</td><td></td><td></td><td></td><td></td><td>19</td><td>2</td><td>14</td><td>5</td><td>26</td><td>7</td><td>11</td><td>5</td><td>34</td><td>15</td><td>24</td><td>15</td></tr><tr><td>2006</td><td></td><td></td><td></td><td></td><td>28</td><td>14</td><td>10</td><td>6</td><td>11</td><td>17</td><td>8</td><td>6</td><td></td><td></td><td></td><td></td></tr><tr><td>2007</td><td>13</td><td>16</td><td>6</td><td>9</td><td>7</td><td>18</td><td>12</td><td>9</td><td>39</td><td>6</td><td>12</td><td>9</td><td>35</td><td>16</td><td>39</td><td>16</td></tr><tr><td>2008</td><td>12</td><td>8</td><td>9</td><td>15</td><td>16</td><td>4</td><td>18</td><td>9</td><td>28</td><td>7</td><td>9</td><td>13</td><td>14</td><td>22</td><td>28</td><td>32</td></tr><tr><td>2009</td><td>34</td><td>10</td><td>13</td><td>15</td><td>29</td><td>11</td><td>19</td><td>12</td><td>30</td><td>7</td><td>14</td><td>12</td><td>10</td><td>6</td><td>39</td><td>27</td></tr><tr><td>2010</td><td>17</td><td>8</td><td>8</td><td>5</td><td>12</td><td>9</td><td>13</td><td>2</td><td>15</td><td>8</td><td>7</td><td>7</td><td>19</td><td>6</td><td>35</td><td>21</td></tr><tr><td>2011</td><td>25</td><td>16</td><td>12</td><td>5</td><td>24</td><td>5</td><td>15</td><td>3</td><td>5</td><td>2</td><td>7</td><td>5</td><td>34</td><td>11</td><td>30</td><td>17</td></tr><tr><td>2012</td><td>16</td><td>21</td><td>18</td><td>10</td><td>31</td><td>14</td><td>13</td><td>5</td><td>12</td><td>10</td><td>5</td><td>4</td><td>32</td><td>17</td><td>25</td><td>11</td></tr><tr><td>2013</td><td>39</td><td>11</td><td>27</td><td>14</td><td>24</td><td>17</td><td>17</td><td>8</td><td>11</td><td>9</td><td>6</td><td>9</td><td>30</td><td>16</td><td>28</td><td>22</td></tr><tr><td>2014</td><td>41</td><td>20</td><td>31</td><td>11</td><td>13</td><td>17</td><td>18</td><td>11</td><td>3</td><td>6</td><td>6</td><td>8</td><td>1</td><td>19</td><td>27</td><td>16</td></tr><tr><td>2015</td><td>50</td><td>23</td><td>26</td><td>9</td><td>25</td><td>7</td><td>19</td><td>11</td><td>16</td><td>2</td><td>7</td><td>6</td><td>2</td><td>1</td><td>23</td><td>15</td></tr><tr><td>2016</td><td>49</td><td>25</td><td>35</td><td>15</td><td>28</td><td>15</td><td>21</td><td>12</td><td>21</td><td>11</td><td>5</td><td>6</td><td>20</td><td>3</td><td>17</td><td>21</td></tr><tr><td>2017</td><td>59</td><td>28</td><td>43</td><td>16</td><td>17</td><td>21</td><td>19</td><td>13</td><td>14</td><td>10</td><td>10</td><td>8</td><td>35</td><td>17</td><td>12</td><td>18</td></tr><tr><td>2018</td><td>45</td><td>31</td><td>45</td><td>21</td><td>25</td><td>13</td><td>28</td><td>10</td><td>14</td><td>9</td><td>14</td><td>8</td><td>12</td><td>24</td><td>20</td><td>22</td></tr><tr><td>2019</td><td>4</td><td>24</td><td>46</td><td>16</td><td>20</td><td>12</td><td>25</td><td>7</td><td></td><td></td><td></td><td></td><td>15</td><td>8</td><td>27</td><td>18</td></tr><tr><td>Total</td><td>625</td><td>399</td><td>767</td><td>368</td><td>631</td><td>318</td><td>531</td><td>227</td><td>312</td><td>128</td><td>139</td><td>125</td><td>340</td><td>203</td><td>467</td><td>293</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>**Figure S1.** Variance of emergence date as a function of the distance between emergence locations, all meadow and individual categories considered. The black dots and line represent the empirical variogram, the size of the point representing the number of pairs of points considered in each distance bin. The red line represents the modelled variogram. The variance stabilized around 40 meters distance between two emergence locations.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure S2.** Correlations between soil temperatures at different depths. Upper part: Daily mean of soil temperatures at various depth (different colors) over the 28 years of study. Lower part: Pearson correlation matrix between soil temperature at depths between surface and two-meter depth (* = p < 0.05, ** = p < 0.01, *** = p < 0.001).</figcaption></figure>
From 2008 to 2019, we were able to acquire air temperatures from the Sheep River weather station (the closest weather station from the study site, with records starting in September 20008, https://acis.alberta.ca/weather-data-viewer.jsp) and compare it to our model estimations at the same location. This was a first validation of the model, in one location not far away from the study site (see figure S3 below), where the microclimate model outputs were highly correlated with empirical data when ran in the place of the weather station. Besides the good predictive ability of the model in other places, this further validated the use of the microclimate model in our study.

1196

<figure><img src="image_23.png" /><figcaption>**Figure S3.** (A) Example of decomposition of a time series of temperature (in °C) recorded at the weather station sheep river between 2008 and 2019. The upper part is the observed signal, and below the trend without the seasonal and random patterns, which can be extracted and analyzed. (B) Correlation of two temperature trends (°C) at sheep river location, between 2008 and 2019 (r = 0.94, Pearson correlation test, p < 0.001). The trends, extracted from decompositions from model estimations (black line) and empirical data (red line), had variations closely related.</figcaption></figure>

1197

1198

1199

1200

1201

1202

1203

1204

1205
<table><caption>**Table S2.** Summary table of spatial variation in emergence date. The mean emergence date (in ordinal day) of each cluster within each meadow (MA, MB, MC, DOT) is presented along with standard errors.</caption><thead><tr><th>Cluster</th><th>Meadow</th><th>Emergence date</th><th>Cluster</th><th>Meadow</th><th>Emergence date</th></tr></thead><tbody><tr><td>bigmiddle</td><td>MA</td><td>120 ± 0.4</td><td>middleeast</td><td>MC</td><td>123 ± 0.9</td></tr><tr><td>bottomeast</td><td>MA</td><td>126 ± 1.2</td><td>middlewest</td><td>MC</td><td>121 ± 0.8</td></tr><tr><td>bottommiddle</td><td>MA</td><td>124 ± 0.7</td><td>overroad</td><td>MC</td><td>122 ± 6.3</td></tr><tr><td>bottomwest</td><td>MA</td><td>122 ± 0.6</td><td>river</td><td>MC</td><td>121 ± 1.4</td></tr><tr><td>extremeast</td><td>MA</td><td>127 ± 3.7</td><td>riverwest</td><td>MC</td><td>120 ± 2.6</td></tr><tr><td>extremewest</td><td>MA</td><td>118 ± 2.6</td><td>roadeast</td><td>MC</td><td>120 ± 0.5</td></tr><tr><td>foresteast</td><td>MA</td><td>124 ± 0.5</td><td>roadwest</td><td>MC</td><td>124 ± 1.1</td></tr><tr><td>middle</td><td>MA</td><td>122 ± 0.6</td><td>bottommiddle</td><td>DOT</td><td>121 ± 0.9</td></tr><tr><td>river</td><td>MA</td><td>119 ± 1.5</td><td>foresteast</td><td>DOT</td><td>123 ± 1.5</td></tr><tr><td>roadeast</td><td>MA</td><td>122 ± 1.1</td><td>forestmiddle</td><td>DOT</td><td>118 ± 1.3</td></tr><tr><td>roadwest</td><td>MA</td><td>120 ± 0.9</td><td>forestwest</td><td>DOT</td><td>121 ± 4.2</td></tr><tr><td>bottomeast</td><td>MB</td><td>118 ± 0.5</td><td>middle</td><td>DOT</td><td>118 ± 0.7</td></tr><tr><td>bottomwest</td><td>MB</td><td>118 ± 0.4</td><td>middleroad</td><td>DOT</td><td>117 ± 0.7</td></tr><tr><td>middleslope</td><td>MB</td><td>121 ± 0.5</td><td>onroad</td><td>DOT</td><td>117 ± 0.9</td></tr><tr><td>smallslope</td><td>MB</td><td>122 ± 2.9</td><td>roadeast</td><td>DOT</td><td>118 ± 1.5</td></tr><tr><td>smallwest</td><td>MB</td><td>130 ± 4.2</td><td>roadeasteast</td><td>DOT</td><td>123 ± 3.2</td></tr><tr><td>top</td><td>MB</td><td>127 ± 1.2</td><td>roadwest</td><td>DOT</td><td>124 ± 1.3</td></tr><tr><td>topeast</td><td>MB</td><td>126 ± 0.9</td><td>smalloverroad</td><td>DOT</td><td>118 ± 5.9</td></tr><tr><td>topwest</td><td>MB</td><td>122 ± 0.7</td><td>topoverroad</td><td>DOT</td><td>114 ± 0.7</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>**Figure S4.** Long-term directional changes in emergence date. Linear temporal trends of annual mean emergence dates (in ordinal day) on the four different meadows (MA, MB, MC, DOT) and for each individual category. Points represent annual means ± standard error in emergence dates of yearlings (green squares), 2-year-old individuals (orange circles), older females (blue triangles), and older males (pink diamonds). Plotted lines correspond to significant temporal trends and their 95% confidence interval.</figcaption></figure>
<table><caption>**Table S3.** Table of the contribution of the 11 climate variables to the principal component analysis performed on annual time series of climate variables per cluster. Loadings of climate factors on the three first components of the principal component analysis (cumulative explained variance: 78.3%).</caption><thead><tr><th>Factor</th><th>Dim.1</th><th>Dim.2</th><th>Dim.3</th></tr></thead><tbody><tr><td>Humidity_inactive</td><td>-0.3</td><td><strong>0.6</strong></td><td>-0.0</td></tr><tr><td>Depthsnow_inactive</td><td>-0.2</td><td><strong>-0.6</strong></td><td>0.5</td></tr><tr><td>Tempair_inactive</td><td>0.3</td><td><strong>0.8</strong></td><td>-0.3</td></tr><tr><td>Wind_inactive</td><td><strong>-0.6</strong></td><td>0.4</td><td><strong>0.6</strong></td></tr><tr><td>Tempsoil_inactive</td><td><strong>0.7</strong></td><td><strong>0.6</strong></td><td>0.2</td></tr><tr><td>Humidity_active</td><td><strong>-0.7</strong></td><td>0.4</td><td>-0.1</td></tr><tr><td>Tempair_active</td><td><strong>0.8</strong></td><td>0.1</td><td>0.5</td></tr><tr><td>Wind_active</td><td><strong>-0.8</strong></td><td>0.4</td><td>0.4</td></tr><tr><td>Tempsoil_active</td><td><strong>0.9</strong></td><td>0.2</td><td>0.3</td></tr><tr><td>Days_snow</td><td><strong>-0.4</strong></td><td>-0.1</td><td>-0.1</td></tr><tr><td>Day_snowmelt</td><td><strong>-0.4</strong></td><td>0.1</td><td>-0.3</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>**Figure S5.** Spatial variability of climate. Principal components, built from climate variables, compared between meadows (highest boxplots of each panel) and clusters of the different meadows (lowest boxplots of each panels).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure S6.** Long-term directional changes in climate. Linear temporal trends of annual mean values of climate variables averaged by meadow on the four different meadows (MA, MB, MC, DOT). Plotted lines correspond to significant temporal trends and their 95% confidence interval (shaded area).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure S7.** Correlations between the different climate variables. The correlation matrix was established from climate variables used in the analyses on the link between climate and emergence date, averaged by year and cluster. Colors correspond to the different clusters (* = p < 0.05, ** = p < 0.01, *** = p < 0.001).</figcaption></figure>
<figure><img src="image_28.png" /></figure>

**Figure S8.** Spatial correlations between the different climate variables. The correlation matrix was established from climate variables used in climate-ED analyses and centered by year (* = p < 0.05, ** = p < 0.01, *** = p < 0.001).
<figure><img src="image_29.png" /><figcaption>**Figure S9.** Temporal correlations between the different climate variables. The correlation matrix was established from climate variables used in climate-ED analyses and centered by cluster (* = p < 0.05, ** = p < 0.01, *** = p < 0.001).</figcaption></figure>