Anouch Tamian
email: anouch.tamian@iphc.cnrs.fr

Vincent A Viblanc

F Stephen Dobson

Peter Neuhaus

Tracey L Hammer

Anna P Nesterova

Shirley Raveh

Amy L Skibiel

David Broussard

Theodore G Manno

Stephen Dobson

Nandini Rajamani

Claire Saraux

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

Keywords: Phenology, climate variability, hibernation, exogenous effects, microclimate model, spatial variability, snow cover

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Walther | Ecological responses to recent climate change[END_REF][START_REF] Martay | Impacts of climate change on national biodiversity population trends[END_REF][START_REF] Pires | Interactive effects of climate change and biodiversity loss on ecosystem functioning[END_REF]. In particular, variations in winter and spring onsets are occurring globally in seasonal ecosystems around the world [START_REF] Menzel | Growing season extended in Europe[END_REF][START_REF] Kunkel | Temporal variations in frost-free season in the United States: 1895-2000[END_REF][START_REF] Kreyling | Vanishing winters in Germany: soil frost dynamics and snow cover trends, and ecological implications[END_REF], Pachauri et al. 2014a, 2014b), affecting the annual timing, known as phenology, of life cycle events (e.g. reproduction) in many animal species [START_REF] Charmantier | Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird Population[END_REF][START_REF] Cohen | A global synthesis of animal phenological responses to climate change[END_REF][START_REF] Sparks | The effect of spring temperature on the appearance dates of British butterflies 1883-1993[END_REF]. 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 andYohe 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 [START_REF] Radchuk | Adaptive responses of animals to climate change are most likely insufficient[END_REF][START_REF] Strode | Implications of climate change for North American wood warblers (Parulidae)[END_REF][START_REF] Visser | Warmer springs disrupt the synchrony of oak and winter moth phenology[END_REF][START_REF] Winder | Climate Change Uncouples Trophic Interactions in an Aquatic Ecosystem[END_REF][START_REF] Chmura | The mechanisms of phenology: the patterns and processes of phenological shifts[END_REF]Parmesan, 2007;Thackeray et al., 2016; see [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Root | Fingerprints of global warming on wild animals and plants[END_REF] for a meta-analyses).

Because global climate change is altering phenological cycles [START_REF] Menzel | Phenology: Its Importance to the Global Change Community[END_REF][START_REF] Cleland | Shifting plant phenology in response to global change[END_REF][START_REF] Chmura | The mechanisms of phenology: the patterns and processes of phenological shifts[END_REF], and because phenological shifts have marked consequences for individual reproduction (e.g. match-mismatch of reproduction and peak resource availability; [START_REF] Durant | Climate and the match or mismatch between predator requirements and resource availability[END_REF][START_REF] Visser | Shifts in phenology due to global climate change: the need for a yardstick[END_REF] and population dynamics [START_REF] Ozgul | Coupled dynamics of body mass and population growth in response to environmental change[END_REF], 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 [START_REF] Graae | On the use of weather data in ecological studies along altitudinal and latitudinal gradients[END_REF][START_REF] Potter | Microclimatic challenges in global change biology[END_REF]. However, individuals are often more affected by local climate conditions, which are shaped by spatial heterogeneity in habitat topography, soil composition, or sun exposure [START_REF] Suggitt | Habitat microclimates drive fine-scale variation in extreme temperatures[END_REF][START_REF] Maclean | Fine-scale climate change: modelling spatial variation in biologically meaningful rates of warming[END_REF][START_REF] Zellweger | Advances in Microclimate Ecology Arising from Remote Sensing[END_REF]. This is especially true in rugged areas, such as mountain habitats, where heterogeneity is high over relatively small spatial scales [START_REF] Huggett | Geoecology: An Evolutionary Approach[END_REF][START_REF] Nagy | The Biology of Alpine Habitats[END_REF][START_REF] Gultepe | Mountain Weather: Observation and Modeling[END_REF]. 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) [START_REF] Barry | Mountain Weather and Climate[END_REF]). Due to a lack of climate data at a fine resolution (<1 km), few studies have examined phenological responses to microclimatic conditions [START_REF] James | The Grinnellian Niche of the Wood Thrush[END_REF][START_REF] Suggitt | Habitat microclimates drive fine-scale variation in extreme temperatures[END_REF][START_REF] Potter | Microclimatic challenges in global change biology[END_REF][START_REF] Varner | The Importance of Biologically Relevant Microclimates in Habitat Suitability Assessments[END_REF][START_REF] Paruchuri | Microhabitat use by plateau pikas: living on the edge[END_REF].

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.) hibernate for ~8 months of the year, and are only active for a short summer period during which reproduction occurs [START_REF] Murie | Annual Variation of Spring Emergence and Breeding in Columbian Ground Squirrels (Spermophilus columbianus)[END_REF], Young 1990b[START_REF] Dobson | Seasonal activity and body mass of Columbian ground squirrels[END_REF]). 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 [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF]. While this study suggested phenological plasticity in ground squirrel ED, likely due to climate change (see also [START_REF] Dobson | Fitness implications of seasonal climate variation in Columbian ground squirrels[END_REF], it only considered adult females from a single population and used widescale 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 [START_REF] Wiggett | Movements of intercolony natal dispersers in the Columbian ground squirrel[END_REF], adult ground squirrels are largely highly philopatric and have a relatively small ranges over their lifetimes (e.g. female range rarely exceeds 1 km 2 and male range around 4.2km²; [START_REF] Arnaud | Philopatry and within-colony movements in Columbian ground squirrels[END_REF][START_REF] Festa-Bianchet | Territoriality in adult female Columbian ground squirrels[END_REF][START_REF] Harris | Inheritance of Nest Sites in Female Columbian Ground Squirrels on JSTOR[END_REF][START_REF] Murie | Territoriality and dominance in male Columbian ground squirrels (Spermophilus columbianus)[END_REF]. 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 andHarris 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 threestep 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 [START_REF] Kearney | How will snow alter exposure of organisms to cold stress under climate warming?[END_REF][START_REF] Kearney | A method for computing hourly, historical, terrain-corrected microclimate anywhere on earth[END_REF][START_REF] Lembrechts | Microclimatic conditions anywhere at any time![END_REF][START_REF] Maclean | Predicting future climate at high spatial and temporal resolution[END_REF][START_REF] Saleeba | Using Biophysical Models to Improve Survey Efficiency for Cryptic Ectotherms[END_REF] 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 [START_REF] Norquay | Hibernation phenology of Myotis lucifugus[END_REF][START_REF] Williams | Phenology of hibernation and reproduction in ground squirrels: integration of environmental cues with endogenous programming[END_REF][START_REF] Fietz | They like it cold, but only in winter: Climatemediated effects on a hibernator[END_REF], there are two theoretical ways that climate could affect ground squirrel emergence dates: 1) by modifying individual energy expenditure during hibernation [START_REF] Humphries | Climate-mediated energetic constraints on the distribution of hibernating mammals[END_REF][START_REF] Pretzlaff | Impact of Climatic Variation on the Hibernation Physiology of Muscardinus avellanarius[END_REF][START_REF] Geiser | Hibernation[END_REF], and 2) through carryover effects from one active season to another [START_REF] Norris | Carry-over effects and habitat quality in migratory populations[END_REF][START_REF] Saino | Sex-dependent carry-over effects on timing of reproduction and fecundity of a migratory bird[END_REF][START_REF] Moore | On the evolution of carry-over effects[END_REF]. 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 andBoag 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 [START_REF] Davis | Hibernation and Circannual Rhythms of Food Consumption in Marmots and Ground Squirrels[END_REF][START_REF] Humphries | Bioenergetic Prediction of Climate Change Impacts on Northern Mammals1[END_REF][START_REF] Fietz | They like it cold, but only in winter: Climatemediated effects on a hibernator[END_REF].

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 [START_REF] Mcguire | Temperature alone is insufficient to understand hibernation energetics[END_REF]. Specifically, we expected areas or years with earlier snowmelt and higher underground temperature to be associated with earlier ED (see [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF]. 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 [START_REF] Shaw | The Storing Habit of the Columbian Ground Squirrel[END_REF].

Overwinter survival depends on fat reserves accumulated during the previous season (e.g. [START_REF] Murie | The Relationship of Body Weight to Overwinter Survival in Columbian Ground Squirrels[END_REF], 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 [START_REF] Dobson | Seasonal activity and body mass of Columbian ground squirrels[END_REF][START_REF] Neuhaus | Timing of Hibernation and Molt in Female Columbian Ground Squirrels[END_REF]. Thus, we expected unfavorable forage conditionse.g. low rainfall (reflected by humidity) to be associated with low primary productivity [START_REF] Dobson | The influence of food resources on life history in Columbian ground squirrels[END_REF][START_REF] Neuhaus | Effects of a late snowstorm and rain on survival and reproductive success in Columbian ground squirrels (Spermophilus columbianus)[END_REF], Dobson et al. 2016), and strong winds to be associated with more time spent by animals in vigilance rather than foraging [START_REF] Fairbanks | Mechanisms of the group-size effect on vigilance in Columbian ground squirrels: dilution versus detection[END_REF]. Through the decrease of fur thermal insulation [START_REF] Mccafferty | Animal thermoregulation: a review of insulation, physiology and behaviour relevant to temperature control in buildings[END_REF], 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 [START_REF] Murie | The Relationship of Body Weight to Overwinter Survival in Columbian Ground Squirrels[END_REF]. 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 [START_REF] Murie | Territoriality and dominance in male Columbian ground squirrels (Spermophilus columbianus)[END_REF][START_REF] Murie | Social interactions and dominance relationships between female and male Columbian ground squirrels[END_REF]. 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 (primarly 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 1527DOT (5.7 ha, from -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. 1,National Band & Tag Co.). Sex was determined based on visual inspection of genitalia [START_REF] Murie | Annual Variation of Spring Emergence and Breeding in Columbian Ground Squirrels (Spermophilus columbianus)[END_REF]. Each year, the meadows were monitored daily from before the first emergences from hibernation (early April each year) to the end of lactation for breeding females. Unfortunately, dates of entry into hibernation were not recorded and we were not able to follow the entire cycle of annual phenology. From mid-April to the end of May each year, ground squirrels were captured within a couple of days of emergence from hibernation (usually on the day of emergence or the following morning) using live-traps (National Live Traps; Tomahawk Co., WI, USA: 13 x 13 x 40 cm³) baited with peanut butter. The first day of observation or capture (usually the same) of a given squirrel was used as a proxy for hibernation emergence date (ED), typically confirmed by the squirrel's appearance and physical condition (large skin flakes from hibernation lost within few days, large abdominal flaps of skin where fat reserves were accumulated and lost, no defecation upon capture; [START_REF] Murie | Annual Variation of Spring Emergence and Breeding in Columbian Ground Squirrels (Spermophilus columbianus)[END_REF]. Emergence was recorded as an ordinal date (i.e. number of days after 1 st January each year). Emergence data were not collected in 2006 on Meadow DOT, in 2003-2006on MA, and in 2002[START_REF] Strode | Implications of climate change for North American wood warblers (Parulidae)[END_REF] on MC. The location of each initial capture was recorded at a ± 1.0m resolution, on 10 x 10 m Cartesian grids of color-flagged wires at each meadow.

Hibernacula are usually used by the animals during a few days after emergence until collective burrows were re-excavated for the active season. First capture locations should thus provide a reasonable approximation of the area where individuals hibernated.

Recorded Cartesian coordinates were converted to GPS coordinates using a database of GPS and Cartesian coordinates collected in 2018-2019. We estimated the error associated with the conversion of Cartesian to GPS coordinates at a median of 5.

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, 2year-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).

Spatial clustering

To evaluate spatial heterogeneity in emergence dates within a meadow, we investigated the spatial autocorrelation of emergence dates. Specifically, we examined whether correlation between emergence dates at two locations depended on the distance between them. To do so, we used variograms [START_REF] Matheron | Principles of geostatistics[END_REF]) that measured the value of the spatial variance at different distance intervals. For any two data points, we measured their spatial distance and estimated the variance in ED. The variance was then averaged over bins of distance intervals. Empirical variograms were based on all emergence locations, but with pairs of points considered only if points were in the same year and within the same meadow.

We modelled variograms using a least square method [START_REF] Rivoirard | Geostatistics for Estimating Fish Abundance[END_REF]. The empirical variogram of emergence date exhibited spatial structure (i.e., increasing variance when the distance between sampling points increased before levelling off; see Appendix S1, Figure S1), meaning that the emergence dates were not random, such that spatially close locations had closer emergence dates than more distant locations. Based on the modelled variogram, we found that variance in emergence date stabilized at a range of about 40 m. Thus, hibernation emergence locations were grouped into spatially relevant clusters according to their GPS coordinates to understand whether individuals on some parts of the meadow emerged earlier than others.

We built clusters by grouping emergence locations according to an iterative algorithm (adapted from [START_REF] Woillez | Indices for capturing spatial patterns and their evolution in time, with application to European hake (Merluccius merluccius) in the Bay of Biscay[END_REF], see also [START_REF] Saraux | Spatial Structure and Distribution of Small Pelagic Fish in the Northwestern Mediterranean Sea[END_REF]. The algorithm started with each location as its own cluster. Each iteration of the algorithm then: (1) calculated the distances between the center of gravity of all clusters on the meadow; (2) merged the two closest clusters if the shortest distance among all cluster distances was below a threshold distance of 40 m (according to above variogram results). The process was reiterated until 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 longterm released version, QGIS.org, QGIS Geographic Information System, 2021) and presented on a Google satellite map © 2021 (see Fig. 1.B).

Microclimate data

We calculated hourly microclimatic conditions using the [START_REF] Kearney | A method for computing hourly, historical, terrain-corrected microclimate anywhere on earth[END_REF] 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, [START_REF] Kearney | A method for computing hourly, historical, terrain-corrected microclimate anywhere on earth[END_REF] downscales global atmospheric climate-forcing data (NCEP Reanalysis, historical data from the National Centre for Environmental Predictions, [START_REF] Kalnay | The NCEP/NCAR 40-Year Reanalysis Project[END_REF], 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; [START_REF] Hollister | elevatr: Access Elevation Data from Various APIs[END_REF], soil composition and hydrological properties (acquired from the 'SoilGrids' database in R; [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF], orientation, and hill shade.

The microclimate model increases the quality and resolution of large-scale weather data, and has been empirically validated [START_REF] Kearney | How will snow alter exposure of organisms to cold stress under climate warming?[END_REF][START_REF] Lembrechts | Microclimatic conditions anywhere at any time![END_REF][START_REF] Maclean | Predicting future climate at high spatial and temporal resolution[END_REF][START_REF] Saleeba | Using Biophysical Models to Improve Survey Efficiency for Cryptic Ectotherms[END_REF]. 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 centimers above ground. Thus, using the microclimate model enabled us to obtain climate data over the entire study period at a heigh 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 modelgenerated 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[START_REF] Neuhaus | Timing of Hibernation and Molt in Female Columbian Ground Squirrels[END_REF], 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.

Data analyses

Quantifying variance in emergence date

Emergence date variance partitioning

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

Emergence date ~ (1|sex) + (1|agecategory) + (1| ID) + (1|meadow/cluster) + (1|year) [1]

This model allowed us to partition variance components between the different random terms;

(1|meadow/cluster) representing both the effect of the meadow, and the effect of the cluster nested within the meadow.

Effects of sex and age on emergence date

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

Emergence date ~ sex + agecategory+ sex*agecategory+ (1|ID) + (1|meadow/cluster) + (1|year)

[2]

Spatio-temporal variation in emergence date

Spatial variation in emergence date

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

Emergence date ~ meadow[MA vs. MB vs. MC vs. DOT] + (1|year) + (1| sex-agecategory) + (1|ID) [3]

Emergence date ~ cluster + (1|year) + (1|sex-agecategory) + (1|ID) [4] (one LMM per meadow)

Temporal variation in emergence date

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

Interannual variation in EDs

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

Emergence date ~ s(year, by = meadow)+ s(year, by = sex-agecategory) + sex-agecategory + meadow + (1|ID) [5]

Long-term directional changes in ED

In a context of climate change, most phenological studies showed advanced spring season [START_REF] Walther | Ecological responses to recent climate change[END_REF][START_REF] Cleland | Shifting plant phenology in response to global change[END_REF][START_REF] Dingemanse | Changing temperature regimes have advanced the phenology of Odonata in the Netherlands[END_REF]). Yet, in Columbian ground squirrels, females older than 3 displayed a delay in emergence over a 20year period [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF]. To test whether this long-term directional pattern was consistent over the 4 populations and over other sex/age categories, we ran the previous model (model [5]) as a LMM, and assessed interactions between continuous (year) and categorical (meadow or sex-agecategory) variables:

Emergence date ~ year + sex-agecategory + meadow + year*meadow + year*sex-agecategory + (1|ID) [6]

Spatio-temporal variation in climate Spatial variation in climate variables

To investigate microclimates, we assessed correlations among the 11 previously defined climate variables (air temperature during the active and inactive period, soil temperature at 1m depth during the active and inactive period, relative humidity during the active and inactive period, wind speed during the active and inactive period, snow depth during inactive period, first day of calendar year without snow and number of days with snow during the active period), using a Principal Component Analysis (PCA) on data averaged by year and by cluster on each meadow. The objective was to identify principal variations related 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 (depthsnow) and soil temperature (tempsoil) at a 1-meter depth during hibernation (inactive phase), and (2) air temperature (tempair), relative humidity and wind speed during the previous summer (active phase), in addition to the ordinal date of spring snowmelt (daysnowmelt) and the number of days with snow during the active period (dayssnow).

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)

Climate change

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

Climateemergence date relationships

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

Mean relationships between climate and emergence date

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

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

Individual responses in emergence date might differ due to different constraints between spatial heterogeneity (related to spatial living locations) and temporal heterogeneity (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:

Emergence date ~ depthsnow_winter + tempsoil_winter + daysnowmelt + dayssnow_active + (1|cluster) + (1|ID) [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.

Emergence date ~ depthsnow_winter + tempsoil_winter + tempair_activen-1 + humidityactiven-1 + windactiven-1 + daysnowmelt + dayssnow_active + (1|year) + (1|ID) [15]

Statistics 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, [START_REF] Box | An Analysis of Transformations[END_REF] to the data (monotonous transformation), or removed extreme values from the data (i.e., emergence dates later than the 22 nd 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 [START_REF] Lenth | Least-Squares Means: The R Package lsmeans[END_REF]. 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, [START_REF] Lê | FactoMineR: an R package for multivariate analysis[END_REF]. 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, [START_REF] Burnham | Multimodel Inference: Understanding AIC and BIC in Model Selection[END_REF]. 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 [START_REF] Zuur | A protocol for data exploration to avoid common statistical problems[END_REF]. 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.

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 (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).

Spatio-temporal variation in emergence date

Spatial variation in emergence date

Analyses on spatial differences in hibernation emergence dates showed that individuals on Meadow DOT emerged 2.2 ± 0.5 days and 1.8 ± 0. S2).

Temporal variation in emergence date

Interannual variability in ED

Emergence date exhibited marked interannual variations over the course of the study on all meadows and for each age-sex category (GAMMs [5]; 2.09 < edf < 8.73, 3.43 < F < 31.61, p < 0.001, Fig. 2), except for older males (≥3-year-old) on Meadow B and C (p = 0.158 and p = 0.652, respectively). Cyclic patterns over time were observed with a period 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)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)

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 (Δ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-yearold individuals (LMMs, ≥3-year-old Females: +0.14 ± 0.

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 variance), was mainly related to wind speed during the winter. The three principal components varied depending on the meadow (LMs [7], p < 0.001), and depending on the cluster (LMs [8], p < 0.001, see Appendix S1, Figure S5).

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

Temporal variation in climate variables

Climate variation

As predicted, all climate variables exhibited temporal variability (GAMs [9] and [10]; 6.51 < edf < 8.78, 2.79< F < 15.33, all p < 0.01, Fig. 4). Interannual variation in these variables was similar among meadows (no year*meadow interactions retained, N = 28 years), although the absolute values in wind speed and relative humidity varied across meadows. This confirmed previous results that Meadow A was the windiest and most humid meadow, while Meadow B was the least humid and windy, Meadows C and DOT being intermediate.

Summer air temperature exhibited particularly high interannual variation, increasing and decreasing from 1 to 4°C from one year to the next (e. g. between 1997 and 1998) 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. observations, N = 1450 individuals, ≥3-year-old females: t = 2.389, p = 0.017, n = 1392 observations, N = 559 individuals). Males of 2-year-old and ≥3-year-old seemed to emerge around 2 days later in locations where the snow depth was 1cm higher during hibernation (LMMs, 2-year-old males: t = 2.303, p = 0.023, n = 280 observations, N = 280 individuals, ≥3-year-old males: t = 2.085, p = 0.039, n = 693 observations, N = 300 individuals, effect only marginally significant after transforming the data, Box-Cox transformation, p = 0.058).

Soil temperature during winter tended to have an effect on the emergence date of yearlings and older males, with individuals emerging earlier from warmer hibernacula (Fig. 6, LMMs, yearlings: t = -1.625, p = 0.105, n = 1450 observations, N = 1450 individuals, ≥3-year-old males: t = -1.498, p = 0.137, n = 693 observations, N = 300 individuals).

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

When conducting the analysis, air temperature during the previous summer had a VIF greater than three (between 4.154 and 5.249) and was removed from the model (probably due to its high positive correlation with soil temperature, see Appendix S1, Figure S9). Regarding temporal effects of climate on emergence date, the ordinal day of snow melt was positively correlated to emergence date for all individuals excepting 2-year-old males (emergence dates observations, N = 300 individuals). Thus, individuals emerged later in years when snow melted later (around 0.1 day later for a day of delayed of snowmelt, see Fig. 7). Similarly, active periods with high wind speed correlated positively with delayed emergence dates of older individuals the next year (between 13 and 17 days delay per m/s increase, LMMs, ≥3year-old females: t = 3.357, p = 0.001, n = 1392 observations, N = 559 individuals; ≥3-yearold 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 onemeter depth).

DISCUSSION

Consistent with previous studies [START_REF] Dobson | Seasonal activity and body mass of Columbian ground squirrels[END_REF][START_REF] Murie | Annual Variation of Spring Emergence and Breeding in Columbian Ground Squirrels (Spermophilus columbianus)[END_REF], ground squirrel age and sex accounted for a substantial amount of variation (29.4%) in ED. Older males emerged first from hibernation, followed by older females, 2-year-old individuals, and finally yearlings. Individual identity had a significant effect in the overall pattern of emergence date. This is consistent with previous knowledge on the heritability of emergence date of Columbian ground squirrels (h² = 0.22 in females and 0.34 in males, [START_REF] Lane | A quantitative genetic analysis of hibernation emergence date in a wild population of Columbian ground squirrels[END_REF].

The magnitude of the individual effect (14.6% of variance) was similar to the temporal variance in ED (16.1% of variance), in spite of the strong individual plasticity in ED from year to year [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF]). These results confirm the importance of including amongindividuals heterogeneity in understanding population dynamics [START_REF] Hamel | General conclusion to the special issue Moving forward on individual heterogeneity[END_REF]. 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 [START_REF] Dobson | The limits of phenotypic plasticity in life histories of Columbian ground squirrels[END_REF][START_REF] Dobson | The influence of food resources on life history in Columbian ground squirrels[END_REF][START_REF] Lane | Phenotypic Plasticity and Local Adaptation in a Wild Hibernator Evaluated through Reciprocal Translocation[END_REF]. 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 [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF][START_REF] Neuhaus | Effects of a late snowstorm and rain on survival and reproductive success in Columbian ground squirrels (Spermophilus columbianus)[END_REF], 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 [START_REF] Dobson | Seasonal activity and body mass of Columbian ground squirrels[END_REF], 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 [START_REF] Lampei | Microclimate predicts frost hardiness of alpine Arabidopsis thaliana populations better than elevation[END_REF][START_REF] Lundblad | Nest microclimate and limits to egg viability explain avian life-history variation across latitudinal gradients[END_REF][START_REF] Wang | Burrow characteristics and ecological significance of Marmota himalayana in the northeastern Qinghai-Tibetan Plateau[END_REF]. Concurrently, spatial microclimate differences were found between clusters, suggesting a potential link between ED and local habitat conditions.

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 [START_REF] Finch | 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[END_REF]but see Ockendon et al., 2013), fish [START_REF] Wilson | Phenological mismatch, carryover effects, and marine survival in a wild steelhead trout Oncorhynchus mykiss population[END_REF], and insects (Dingemanse andKalkman 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 [START_REF] Harrison | Carry-over effects as drivers of fitness differences in animals[END_REF]. Climate directly affects food availability/quality or/and foraging patterns [START_REF] Levy | Foraging Activity Pattern Is Shaped by Water Loss Rates in a Diurnal Desert Rodent[END_REF], 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 [START_REF] Fairbanks | Mechanisms of the group-size effect on vigilance in Columbian ground squirrels: dilution versus detection[END_REF], 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 [START_REF] Dobson | The influence of food resources on life history in Columbian ground squirrels[END_REF][START_REF] Neuhaus | Effects of a late snowstorm and rain on survival and reproductive success in Columbian ground squirrels (Spermophilus columbianus)[END_REF], 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 [START_REF] Fietz | They like it cold, but only in winter: Climatemediated effects on a hibernator[END_REF].

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 emerge. In addition, a melting and wet snow pack makes it hard for the animals to easily move around the meadow (VAV, FSD, PN; personal observations) and likely increases the costs of thermoregulation, increasing overall energy expenditure at a time where their remaining fat stores are likely depleted. Local climate may thus act as a strong cue on the timing of ground squirrel emergence, though the question of the immediate mechanism that stimulates emergence remains for future research.

We found a positive spatial correlation between snow depth and emergence of 2-yearold and ≥3-year-old males, suggesting that males of 2-year-old and older emerged later in locations where snow depth was higher during the inactive period. Higher snowpack or delayed snowmelt during the inactive period may also insulate and stabilize soil temperatures and hibernation conditions, decreasing fat reserve consumption, and thus allowing a longer hibernation or emergence in better condition (the latter important for reproduction; [START_REF] Rubach | Testing the reproductive and somatic trade-off in female Columbian ground squirrels[END_REF]. Snowpack thickness is an important variable influencing hibernating species, such as alpine marmots [START_REF] Canale | Differential plasticity of size and mass to environmental change in a hibernating mammal[END_REF], hoary marmots [START_REF] Patil | Winter weather versus group thermoregulation: what determines survival in hibernating mammals?[END_REF], or northern Idaho ground squirrels [START_REF] Goldberg | Hibernation behavior of a federally threatened ground squirrel: climate change and habitat selection implications[END_REF]. A negative correlation between the soil temperature during inactive period (at one-meter depth) and the emergence date of yearlings and ≥3-year-old males also seemed to apply. As for northern Idaho ground squirrels [START_REF] Goldberg | Hibernation behavior of a federally threatened ground squirrel: climate change and habitat selection implications[END_REF], this result possibly indicates that unfavorable hibernating conditions lead to earlier emergence. For example, elevated soil temperatures prevent animals from lowering their body temperature sufficiently for efficient energy conservation during hibernation. In short, fat-storing and hibernating species might be sensitive to local environmental conditions and especially elevated temperatures, which are directly influencing energy expenditure during this period of fast.

Climatic variables correlating with ED differed between sex-age categories.

Explanations for such diverging responses may be both ecological and physiological. Young (1990a) reported the existence of long post-hibernation euthermic phases only in adult males, remaining several days underground before emerging despite favorable aboveground conditions, possibly to allow for gonadal maturation [START_REF] Williams | Sex-Dependent Phenological Plasticity in an Arctic Hibernator[END_REF]. Physiological and phenological differences in sex and age categories might therefore reflect distinct constraints between reproductively mature individuals and non-mature individuals. Yearlings and 2-yearolds that emerge later in the spring likely do so when resources become more abundant and perhaps to avoid predator exposure or the social aggression related to reproduction (Murie andHarris 1982, Constant et al. 2020). By contrast, mature individuals may optimize their emergence date according to the trade-off between energy resources (remaining fat stores, food availability) and the advantages of an early reproduction [START_REF] Neuhaus | Timing of Hibernation and Molt in Female Columbian Ground Squirrels[END_REF]. Furthermore, reproductive strategies might differ between males and females, with males taking advantage of earlier hibernation emergence to seek more breeding opportunities (Richardson's ground squirrels, [START_REF] Michener | Spring emergence schedules and vernal behavior of Richardson's ground squirrels: why do males emerge from hibernation before females?[END_REF]) and establish their territories [START_REF] Manno | Why are Male Columbian Ground Squirrels Territorial?[END_REF]. Thus, environmental influence on a population might differ importantly between individuals with different life history constraints (e.g. age classes; [START_REF] Bonamour | Agedependent phenological plasticity in a wild bird[END_REF][START_REF] Pardo | Evidence for an agedependent influence of environmental variations on a long-lived seabird's life-history traits[END_REF].

Long-term changes in climate in relationship to emergence date

Previous studies have reported shifts in animal phenology linked to changes in climate, species especially responding to temperature (more so at high latitudes) and precipitations [START_REF] Cohen | A global synthesis of animal phenological responses to climate change[END_REF]). In the context of climate change, advancements in spring phenology have been found in the vast majority of cases [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF]; but see [START_REF] Radchuk | Adaptive responses of animals to climate change are most likely insufficient[END_REF]. However, in studying long-term effects of climate change on animal phenology, mammals are poorly represented [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF]but see Radchuck et al. 2019). Among those studied, Marmotini species appear to respond differently to directional changes in climate. The phenology of Uinta ground squirrels did not change over time, despite a locally changing climate [START_REF] Falvo | Seasonal climate effects on the survival of a hibernating mammal[END_REF]). Yellow-bellied marmot have been found to emerge earlier in response to warming air temperatures [START_REF] Inouye | Climate change is affecting altitudinal migrants and hibernating species[END_REF], while Columbian ground squirrels had been shown to display a trend toward later emergence, related to late spring snowstorms [START_REF] Lane | Delayed phenology and reduced fitness associated with climate change in a wild hibernator[END_REF].

Interestingly, our results are only partially consistent with previous findings on Columbian ground squirrels. We found evidence for delayed snowmelt ( 13 Whether variation in ground squirrel phenology is linked to behavioral plasticity versus longterm 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 [START_REF] Lane | Phenotypic Plasticity and Local Adaptation in a Wild Hibernator Evaluated through Reciprocal Translocation[END_REF].

Many species responses to climate are likely to occur at micro-scale, rather than mesoor macro-scale (e.g. sockeye salmon; Onchorhynchus nerka; [START_REF] Martins | Climate effects on growth, phenology, and survival of sockeye salmon (Oncorhynchus nerka): a synthesis of the current state of knowledge and future research directions[END_REF]. 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 climatephenology 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.

  5 ± (SD) 3.4m (n = 225 coordinates). Coordinates were collected for 74% of the emergences, 4351/5873. The missing emergence sites were before 1994 on Meadow A, 1999 on Meadow B, 2006 on Meadow C, and 2012 on Meadow DOT. Emergence locations were also not recorded in 2000-2007 on MA.Thus, sample sizes varied among years.

GAMs [ 9 ]

 9 & [10] as LMs [11] & [12], with each climate variable as a dependent variable and the year, meadow (or cluster) and their interaction as independent variables: Climate variable ~ year*meadow + meadow [11] Climate variable ~ year*cluster + cluster [12] (one LM per meadow)

  4 days earlier than individuals on Meadow C and Meadow B, respectively. Meadow C and Meadow B individuals emerged 3.3 ± 0.4 days and 3.6 ± 0.3 days earlier than Meadow A individuals (LMM [3], post-hoc LSM: all p < 0.001, n = 5873, N = 2436 individuals, 28 years). Similarly, within-meadow differences in emergence dates were found between different clusters (see cluster colors on Fig. 1.B). As an example, on Meadow B, individuals in "bottom-west" cluster emerged between 1.9 ± 0.5 and 8.8 ± 2.5 days earlier than individuals from other clusters (LMM [4], post-hoc LSM: p < 0.001, n = 1315, N = 579 individuals, 21 years). In contrast, individuals from the "top" cluster emerged between 4.6 ± 0.9 and 8.1 ± 0.8 days after individuals from other clusters (post-hoc LSM: p < 0.001, n = 1315, N = 579 individuals, 21 years) (see Fig. 1.B and Appendix S1, Table

  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.

  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).

  , fluctuating on average between 7-13°C. During 2006-2007, air temperature was particularly high, and much warmer than the next four years. Conversely, the relative humidity during summer was 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.

  and climate variables centered by cluster, LMMs [15], yearlings: t = 3.015, p = 0.003 , n = 1450 observations, N = 1450 individuals; 2-year-old females: t = 2.821, p = 0.006, n = 525 observations, N = 525 individuals; 2-year-old males: t = 1.275, p = 0.211, n = 280 observations, N = 280 individuals; ≥3-year-old females: t = 3.246, p = 0.002, n = 1392 observations, N = 559 individuals; ≥3-year-old males: t = 3.610, p < 0.001, n = 693

  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).

Figures

  Figures

Fig. 1 .Fig. 2 .

 12 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.

Fig. 3 .Fig. 4 .

 34 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.

Ordinal

Fig. 5 .

 5 Fig. 5. Climate effects on hibernation emergence date of different individual categories.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).

Fig. 6 .Fig. 7 .

 67 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.

Figure S2 .

 S2 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).

Figure S5 .

 S5 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).

Figure S6 .

 S6 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).

Figure S7 .

 S7 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).

Figure S8 .

 S8 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).

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. 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.

Climateemergence date relationships Mean relationships between climate and emergence date

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

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

Once centered per year, no strong correlations were highlighted between climate variables, and we included all four of them in the model (VIFs < 1.308, see Appendix S1, Figure S8). Day of snowmelt was related positively to emergence date, but only in yearlings and older (≥3-year-old) females; individuals emerged later when they emerged from locations where snow melted later (Fig. 6; ED around 0.2 day later for a day of delayed snowmelt, ED 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). 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.

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.

Table S3. Table of the