Mathilde Poirier
email: mathildepoirier05@gmail.com

Gilles Gauthier

Florent Domine

Dominique Fauteux

Lemming winter habitat: the quest for warm and soft snow

Keywords: Climate change, habitat use, rodent, small mammals, snowpack

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

Introduction

The arrival of below-zero temperatures and of a snow-covered landscape in autumn marks the beginning of a challenging period for boreal and arctic species. Snow, depending on its properties, can either be an ally or a foe for animals living in these environments. For instance, willow ptarmigans (Lagopus lagopus) can benefit from deep snow accumulation as it allows them to browse buds higher up on shrubs and increase their food intake [START_REF] St-Georges | Winter habitat use by ptarmigan, snowshoe hares, red foxes, and river otters in the boreal forest -tundra transition zone of western Quebec[END_REF].

For other herbivores such as reindeer (Rangifer tarandus) and muskoxen (Ovibos moschatus), hard snow or ice crusts resulting from rain-on-snow events can reduce access to ground vegetation and lead to population-wide reduction in survival [START_REF] Rennert | Soil thermal and ecological impacts of rain on snow events in the circumpolar arctic[END_REF][START_REF] Hansen | Climate, icing, and wild arctic reindeer: past relationships and future prospects[END_REF]. Snow properties vary from year to year but also within the landscape. Depending on their specific needs and characteristics, each species should theoretically use the habitat in a way that maximizes snow attributes favoring them. For instance, mammals like red foxes (Vulpes vulpes) or mule deer (Odocoileus hemionus) minimize travelling in deep and soft snow, probably to reduce their energy expenditure for locomotion [START_REF] Halpin | Influence of snow depth on prey availability and habitat use by red fox[END_REF]Bissonette 1988, Coe et al. 2018).

Lemmings are small arctic rodents that benefit from the snowpack for thermoregulation and protection against predators (Duchesne et al. 2011a). These small mammals are the main prey of many predators and these interactions likely lead to their cyclic population fluctuation [START_REF] Gilg | Cyclic dynamics in a simple vertebrate predator-prey community[END_REF][START_REF] Fauteux | Top-down limitation of lemmings revealed by experimental reduction of predators[END_REF]. In winter, lemmings stay active in the snowpack where they dig a network of tunnels to find plants on which they feed or to disperse. They build nests for warmth where they can also reproduce if they have enough energy [START_REF] Millar | On reproduction in lemmings[END_REF], Duchesne et al. 2011a).

Lemmings typically build their nest within a deep snowpack (60 cm or more) where subnivean temperatures reach an optimum (Duchesne et al. 2011a[START_REF] Reid | Lemming winter habitat choice : a snow-fencing experiment[END_REF]. However, in a collared lemming (Dicrostonyx groenlandicus) population in Greenland, use of shallower snowpacks increased as the population density increased, suggesting density-dependent habitat selection resulting from resource competition [START_REF] Schmidt | Long -term patterns in winter habitat selection, breeding and predation in a densityfluctuating, high Arctic lemming population[END_REF]. Most studies have focused on snow depth to explain lemming distribution in the landscape, but other snow properties could also influence their habitat use in winter. For instance, it has been shown that lemmings predominantly use the soft, low-density basal snow layer to move through the snowpack [START_REF] Poirier | What guides lemmings movements through the snowpack?[END_REF], which can vary at the landscape level.

The High Arctic snowpack is shallow, typically 25 cm thick, but often as little as 15 cm (Sturm andBenson 2004, Domine et al. 2012) and is generally comprised of two main layers: a basal soft, fragile layer of low density snow topped by a hard dense wind slab (see Figure 1 of Domine et al. 2018a). It should be noted that this description excludes arctic regions strongly affected by oceanic currents. The basal layer mostly forms at the beginning of winter, when the high temperature gradient between the ground which is still warm and the cold polar air generates strong upward water vapor fluxes that lead to the formation of large, loosely bonded cup-shaped crystals called depth hoar [START_REF] Marbouty | An experimental study of temperature-gradient metamorphism[END_REF]. The decrease in temperature gradient throughout the winter subsequently limits further depth hoar formation and snow evolution is instead governed by wind effects. A greater depth hoar fraction in the snowpack is therefore favored by greater soil moisture which delays freezing, a shallow snowpack, and low wind exposure (Domine et al. 2018a). In areas with a heterogeneous topography (e.g., hummocks, river banks), snow accumulation is greater due to wind redistribution [START_REF] Pomeroy | Physical properties of snow[END_REF]. The presence of shrubs can also increase snow accumulation and depth hoar fraction [START_REF] Sturm | Snow -shrub interactions in arctic tundra : A hypothesis with climatic implications[END_REF]) because branches trap wind-blown snow and limit erosion by strong winds [START_REF] Domine | The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime[END_REF].

Several variables characterize the quality of the different snow layers, as relevant to lemming habitat. Snow density (ρ) is the mass/volume ratio of a snow layer [START_REF] Conger | Instruments and Methods: Comparison of density cutters for snow profile observations[END_REF].

Snow thermal conductivity (keff) is inversely proportional to the thermal insulation properties of the snow and is positively correlated to both density and hardness [START_REF] Domine | Linking the effective thermal conductivity of snow to its shear strength and density[END_REF]. Snow specific surface area (SSA) is a measure of the surface/volume or surface/mass ratio of sampled snow grains and is inversely related to grain size [START_REF] Gallet | Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm[END_REF]). The basal depth hoar is usually relatively soft (ρ = 130-250 kg m -3 , keff = 0.025-0.1 W m -1 K -1 , SSA = 9-11 m 2 kg -1 ; Sturm andBenson 1997, Domine et al. 2018). However, when temperature increases above 0 °C, partial melting and subsequent refreezing of the snow, or a rain-on-snow event, can occur and create hard melt-freeze layers (ρ = 364-480 kg m -3 , keff = 0.11-0.39 W m -1 K -1 , SSA = ~3 m 2 kg -1 ; Sturm et al. 1997, Domine et al. 2009). Although such layers can subsequently metamorphize into depth hoar, it will be fairly hard. In regions where the climate is wetter and milder, meltfreeze layers are more frequent, and warmer air temperature and deeper snowpack will slow down the depth hoar metamorphism [START_REF] Marbouty | An experimental study of temperature-gradient metamorphism[END_REF][START_REF] Sturm | A seasonal snow cover classification system for local to global applications[END_REF]. Another mode of formation of hard depth hoar is when a wind slab forms early in the season. A wind slab is usually dense and hard (ρ = 350-488 kg m -3 , keff = 0.16-0.45 W m -1 K -1 , SSA = 20-30 m 2 kg -1 ; Sturm et al. 1997;[START_REF] Domine | The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime[END_REF]) but when formed early in the season and subjected to a high temperature gradient, it will also metamorphize into hard depth hoar. Hard depth hoar, formed from either refrozen layers or a wind slab, is called indurated depth hoar (ρ = 250 to >350 kg m -3 , keff = 0.05 to > 0.30 W m -1 K -1 ; Domine et al. 2018). Depth hoar with a SSA in the lower range (~8 m 2 kg -1 ) usually indicates formation from a melt-freeze layer while a SSA value ~12 m 2 kg -1 usually indicates formation from a wind slab layer [START_REF] Domine | Three examples where the specific surface area of snow increased over time[END_REF].

With climate change, a delayed onset of the snowpack is expected [START_REF] Liston | The changing cryosphere: Pan-Arctic snow trends (1979-2009)[END_REF] and may reduce the period when lemmings can benefit from it [START_REF] Gilg | Climate change and cyclic predator-prey population dynamics in the high Arctic[END_REF]. Global warming also threatens to increase the frequency of rain-on-snow events, leading to a hardening of the snowpack (Liston andHiemstra 2011, Hansen et al. 2014). As lemmings have greater difficulty digging in harder snow [START_REF] Poirier | Snow hardness impacts intranivean locomotion of arctic small mammals[END_REF]), this could have population-wide impacts by increasing their energy expenditure and reducing their chances of reproduction. Some studies already reported negative consequences of rain-on-snow on lemming and vole populations [START_REF] Ims | Collapsing population cycles[END_REF][START_REF] Kausrud | Linking climate change to lemming cycles[END_REF], Domine et al. 2018b), but a better understanding of their winter habitat use is necessary to assess their vulnerability to future snow conditions.

In this study, we examined winter habitat use and the impact of habitat and snow properties on the probability of reproduction of the brown (Lemmus trimucronatus) and northern collared lemming in the High Arctic. Habitat use was determined by comparing density of lemming winter nests among four habitats (mesic, riparian, shrubland and wetland) and their reproductive activity by comparing the proportion of reproduction in those nests over 13 years. We hypothesized that habitat use and reproduction rates of lemmings would be higher in habitats where snow conditions are expected to facilitate their thermoregulation and locomotion. We predicted that winter nest density and proportion of nests with reproduction would be higher (1) in habitats with a deep snowpack and early snow accumulation because these conditions should offer good thermal insulation, and (2) where the basal snow layer is characterized by a soft and well-developed depth hoar layer because it facilitates locomotion through the snow. We also hypothesized that winter nest density would increase in less favorable habitats in years of high lemming abundance due to a density-dependent spillover effect. Finally, we hypothesized that brown and collared lemmings would use winter habitats differently because collared lemmings are more efficient at digging into the snow than brown lemmings [START_REF] Poirier | Snow hardness impacts intranivean locomotion of arctic small mammals[END_REF].

Material and methods

Study area

This study took place in the Qarlikturvik valley of Bylot Island, Nunavut (73°08'N, 80°00'W), from 2007 to 2019 (Supplementary Material S1 Fig. S1). Winter, which is defined here as the period with snow cover, typically starts in early October and ends in early June [START_REF] Domine | Meteorological, snow and soil data (2013-2019) from a herb tundra permafrost site at Bylot Island, Canadian high Arctic, for driving and testing snow and land surface models[END_REF]. The coldest month is February with a mean temperature of -36.7 °C [START_REF] Domine | Meteorological, snow and soil data (2013-2019) from a herb tundra permafrost site at Bylot Island, Canadian high Arctic, for driving and testing snow and land surface models[END_REF].

Climate in this area is typical of cold and dry Arctic regions unaffected by oceanic currents. For this study, we separated the study area (~51 km 2 ) in four main habitats for lemmings in winter, namely, mesic (67 %), wetland (17 %), riparian (10 %) and shrubland (6 %) (Supplementary Material S1 Fig. S1-S2). The mesic habitat is mostly found upland, in sites with moderate to good drainage. It is characterized by hummocky tundra and the main plant species encountered are Salix arctica, Cassiope tetragona, graminoids and mosses [START_REF] Audet | Feeding ecology of greater snow goose goslings in mesic tundra on Bylot Island, Nunavut, Canada[END_REF]). Wetland habitat (wet) is found in low-lying areas close to sea-level and is characterized by shallow ponds and polygon tundra created by the growth of ice wedges in the permafrost. The vegetation mostly consists of graminoids, such as Carex aquatilis, Eriophorum sheuchzeri and Dupontia fisheri growing through an extensive moss cover [START_REF] Gauthier | The tundra food web of Bylot Island in a changing climate and the role of exchanges between ecosystems[END_REF]. Riparian habitat is restricted to the area along streams or gullies with steep slopes on either side (i.e., 10-30 m from them), running mostly through the mesic habitat or along hills, which allow the formation of snowdrifts in winter [START_REF] Pomeroy | Physical properties of snow[END_REF]. This habitat has a similar vegetation composition to mesic, but near the center of gullies, vegetation becomes more similar to wetland due to higher ground humidity. Finally, shrubland habitat (shrub) resembles mesic but is the only habitat with erect vegetation, primarily Salix richardsonii, rising 10 to 40 cm above the ground [START_REF] Domine | The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime[END_REF]). This habitat is very limited spatially and mostly occurs in the most inland part of the valley, at the base of alluvial fans which provide abundant water supply by channeling snowmelt water.

Lemming demographic parameters

Lemming winter nests were sampled along 20 permanent transects of 500 m each in mesic, riparian and wet habitats between 2007 and 2019 (total of 60 transects each year; Supplementary Material S1 Fig S1). For the shrub habitat, the sampling took place between 2017 and 2019 along 25 transects of 200 m (transects were shorter because that habitat occurred in small-size patches).

These transects were randomly positioned in the landscape from predetermined coordinates and had a predominantly N-S orientation to avoid transects crossing each other. Most of them were spaced by at least 200 m, which ensured spatial independence considering the small home range of lemmings (Banks et al. 1975, Predavec andKrebs 2000). We walked all transects soon after complete snow melt and noted every winter nest detected from the transect line. We measured the perpendicular distance between each nest and the transect and destroyed them so only new nests could be sampled the next year. We identified the lemming species that occupied every nest found based on the size, shape and color of feces [START_REF] Maclean | Cycles in arctic lemmings: winter reproduction and predation by weasels[END_REF][START_REF] Soininen | Highly overlapping winter diet in two sympatric lemming species revealed by DNA metabarcoding[END_REF].

Feces were easy to identify because we only sampled fresh nests from the year. In some cases (about 5 %), feces from both brown and collared lemmings were found, suggesting a mixed occupation. In those cases, we duplicated the nest in our dataset so it would be considered twice in the analysis, once for each species.

We estimated nest density of each species separately in each habitat and year with the distance sampling method [START_REF] Miller | Distance sampling in R[END_REF]. Data from all transects were pooled within each habitat to obtain a single estimate of nest density per habitat annually. We assumed that detection probabilities decreased with distance separating the nest from the transect. For each habitat and year, we modeled the detection probability with different probability distribution functions (half normal, hazard rate, uniform; [START_REF] Buckland | Advanced Distance Sampling: Estimating Abundance of Biological Populations[END_REF]. It was also possible to specify an adjustment factor (i.e., cosine, Hermite polynomial or simple polynomial) to improve the fit of the model. We used the second-order Akaike criterion (AICc) to determine which probability distribution function provided the best fit for each dataset. The model estimated the effective area surveyed by the observer as well as the number of undetected nests based on the detection probability function. Nest density was obtained by the sum of detected and undetected nests, divided by the effective area surveyed in each habitat. Those analyses were performed with the Distance package version 1.0.3. [START_REF] Miller | Distance sampling in R[END_REF] implemented in the R software version 4.1.0 (R Core Team 2021). When models could not be run or performed poorly due to low sample size (about 25% of the time), we obtained a nest density by dividing the number of nests found in the habitat by the total area surveyed assuming a perfect detection 5 m on either side of the transect.

For each lemming nest, we determined the occurrence of reproduction based on the presence of a high number of small-size feces (i.e., at least one third of all feces), which indicates that juveniles once occupied the nest (Duchesne et al. 2011b). In cases of mixed nest occupation, the layering order of feces in the nest indicated which species had reproduced. We estimated the proportion of nests with reproduction for each species, year and habitat by dividing the number of nests with presence of reproduction by the total number of nests found. In years of low lemming densities, the number of winter nests found along transects was very low, which reduced the precision of the proportions of nests with reproduction. In those years, we increased sample size by including winter nests found opportunistically in the field. Opportunistic nests were found while observers conducted other field activities across the study area and were considered randomly sampled. The habitat where opportunistic nests were found was assigned in the field and they were analyzed in the same way as those found along transects. However, they were not used in the density estimation.

Snow physical properties

We sampled snow physical properties in the four habitats by digging snow pits in May before snow melt from 2014 to 2019 except in 2016 due to logistical constraints. We dug between 1 to 11 snow pits in every habitat each year. It was not possible to perform snow pits at the exact same location than the winter nests transects, but we selected sites with similar habitat characteristics within the area covered by our transects (Supplementary Material S1 Fig. S1). In every snow pit, we sampled physical properties of the basal layer of snow between 0 to 5 cm from the ground level. We measured snow density by weighing a fixed volume of snow with a 100 cm 3 box-cutter [START_REF] Conger | Instruments and Methods: Comparison of density cutters for snow profile observations[END_REF]. As a proxy for snow hardness, we measured the thermal conductivity of snow (keff) with a TP02 heated needle probe (Hukseflux Thermal Sensors, Delft, The Netherlands). A needle was inserted in the basal snow layer, avoiding touching the ground, heated with a constant power during 100 s and the temperature at the center of the needle was recorded. The relationship between temperature at the center of the needle and time (on a log-scale) depends on heat dissipation in the environment and can be used to calculate thermal conductivity (Morin et al. 2010[START_REF] Domine | Linking the effective thermal conductivity of snow to its shear strength and density[END_REF]. We used the DUFISSS instrument (Dual Frequency Integrating Sphere for Snow SSA measurement) to measure the snow specific surface area (SSA), which is the surface area per unit of mass [START_REF] Gallet | Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm[END_REF]). The idea consists in illuminating a given volume of snow with a laser diode at 1310 nm to measure the reflected light with a photodiode, and to convert the reflectance measurement into SSA (more details in [START_REF] Gallet | Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm[END_REF][START_REF] Domine | Structure, specific surface area and thermal conductivity of the snowpack around Barrow, Alaska[END_REF]. Finally, we visually determined the depth hoar layer in the snowpack, and we measured the maximal height of this layer and the total height of the snow pit.

We measured mean snow depth in each habitat in mid-May, which is usually the period of maximum snow depth. In mesic, shrub and wet habitats, this was measured using a snow probe along 4 random transects positioned in a square shape, totaling 100 to 200 points per habitat spaced out by 2 m. For riparian habitat, because the snow depth varies considerably at the meter scale within snow drifts, we used the maximum depth recorded in May at three permanent poles centrally located in selected snowdrifts.

From 2016 to 2021, three automated stations (one per habitat) monitored snow temperature with thermistors in riparian, shrub and wet habitats throughout the winter. Since the snow depth is similar between wet and mesic habitats, we assumed that their thermal profiles are also relatively similar even though we recognize that differences in slope, aspect, soil properties and vegetation may affect them. The thermistors were installed at different heights in the snowpack between 0 and 35 cm. We used the temperature at 2 cm above the ground, obtained either from direct measurement at that height or from linear interpolation using the nearest 2 thermistors. We obtained the daily temperature as well as the daily temperature fluctuation (i.e., maximalminimal daily temperature) of the basal snow layer. Due to malfunction of the equipment in some years, we could only compare the three habitats over three winters (2016-17, 2019-20 and 2020-21). Air temperature was measured at two automated stations (shrub, wet) with a ventilated sensor at 2.3 m height [START_REF] Domine | Meteorological, snow and soil data (2013-2019) from a herb tundra permafrost site at Bylot Island, Canadian high Arctic, for driving and testing snow and land surface models[END_REF].

We obtained dates of permanent snow onset from 2017 to 2021 using automatic cameras (Reconyx ® ) that recorded daily photographs of a representative area of each habitat. One camera positioned about 10 m above the ground took pictures in mesic, riparian and wet habitats simultaneously, and a second camera recorded pictures in the shrub habitat. Snow onset was defined as the first day in autumn when snow cover reached or exceeded 80 % of the area and did not return below 50 %. Note that throughout the paper, the year of a given winter is referred to by the year of the spring when most measurements were taken (e.g., winter 2016-2017 is referred to as winter 2017).

Statistical analyses

We analyzed lemming winter nest density with linear mixed-effects models to assess the influence of either habitat or snow parameters, lemming species and their interactions when relevant. The sampling unit used in this analysis was annual nest density estimated in each habitat with the distance sampling method (see Lemming demographic parameters above). We log-transformed (natural log) nest densities to meet normality and homoscedasticity assumptions.

The temporal autocorrelation in model residuals caused by the cyclic dynamics of these rodents was removed by adding year as a random factor. To test for a density-dependent effect on habitat use, we repeated the global model with data divided into years of low and high nest density and compared the effect size and significance of the habitat covariates. The threshold used to separate low and high density was the 2.9 nests/ha, which corresponds to the median of yearly nest density for brown and collared lemmings combined (Supplementary Material S1 Fig. S4).

For datasets including the shrub habitat (available only for 2017-2019), we could not test for a density-dependent effect due to the short time series. To test the influence of snow physical properties on nest densities (2014-2019), we used different combinations of additive parameters (i.e., density and thermal conductivity of the basal snow layer, maximal height of depth hoar, snow depth) and we evaluated the strength of support of each model using the second-order Akaike criterion (AICc). We used model-averaging when several had reasonable statistical support (ΔAICc < 4).

We analyzed the proportion of winter nests with reproduction with binomial models to examine the influence of habitat and species. For the longest time series (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), we used generalized linear mixed quasi-binomial models to overcome overdispersion in the data (MASS package version 7.3-57.1; Ripley et al. 2022), with year included as random factor. For the dataset including shrub (2017-2019), we used binomial models without year as a random effect due to the small dataset. Similarly, we used binomial models to test for influence of species and snow parameters on proportion of nests with reproduction (2014-2019). We evaluated the strength of support of each model using ΔAICc and determined significance of relationships based on the 95% confidence interval of slope parameters.

To establish potential differences in snow properties between the four habitats, we used linear mixed effects models with habitat as the fixed effect and year as a random effect (2014-2019).

Response variables were density, thermal conductivity and specific surface area (SSA) of the basal snow layer, snow depth and maximal height of the depth hoar. We also examined the differences in temperature (with and without correcting for differences in air temperature) and daily temperature fluctuation of the basal snow layer between habitats with data from the automated stations. To avoid autocorrelation problems of repeated measures, we used mean monthly temperatures and mean of daily fluctuations in temperature calculated over 15-day intervals (time intervals were determined to minimize autocorrelation in the data and to maximize sample size). For the snow onset date, we examined the difference only between the wet and riparian habitats because dates were the same for mesic and wet habitats and almost the same for riparian and shrub habitats. In some models we had to perform a natural log, inverse or square-root transformation of the response variable to improve normality and homoscedasticity.

Due to the different scales of variables, we could not use the same transformation consistently.

As a complementary analysis, we explored the degree of association amongst snow variables through a principal component analysis (PCA) and examined variations in PC scores between habitats (see details in Supplementary Material S2).

When relevant, the proportion of variation explained by the models were calculated with the MuMIn package version 1.47.1 (Barton 2022) following Nakagawa and Schielzeth 2013, with R 2 m being the amount of variation explained by fixed factors and R 2 c by both fixed and random factors. In tables, slope parameters (β) are presented with their 95 % CI and differences between groups are calculated with the Tukey's Honestly Significant Difference test (HSD). When logtransformed data where used in the models, we obtained the mean estimates on the linear scale by applying an exponential transformation of the sum of the estimate and half of the variance [START_REF] Feng | Log-transformation and its implications for data analysis[END_REF]) and we calculated the approximate 95% CI with the Cox's method [START_REF] Chami | On efficient confidence intervals for the log-normal mean[END_REF].

Results

Winter nest density

Nest density of both lemming species varied considerably among years, as expected based on their known cyclic population fluctuations, and between habitats (Supplementary Material S1 Fig. S3, Fig. S4). We found evidence that nest densities were highest in the riparian habitat, lowest in the wet habitat and intermediate in the mesic habitat (Fig. 1, Table 1a). Nest densities of collared lemmings were also 2.6 times lower than that of brown lemmings (β = -0.96, CI = [-1.38, -0.54]) and a model with an interaction between lemming species and habitat had less support (ΔAICc = 1.2). We did not find evidence of density-dependent effects in habitat use since models where years were separated in low and high density yielded similar results (i.e., similar slopes) to the global model (Supplementary Material S1 Fig. S5, Table S1). Results of the analysis using years with shrub habitat data were very consistent with the previous analysis despite the smaller dataset (Table 1a). Nest density in the shrub habitat was 1.2 and 3.7 times higher than in the mesic and wet habitats, respectively but 3.8 times lower than in the riparian habitat (Fig. 1 c-d, Table 1a).

When examining the influence of snow parameters on nest density, we found support for a positive effect of snow depth and maximal height of depth hoar on nest density and slightly lower nest densities of collared than brown lemmings (Table 2a, Fig. 2, Supplementary Material S1 Table S2). A model with an inverse transformation of snow depth had less support, providing little evidence of a threshold effect of snow depth on lemming nest density.

Winter reproduction

The proportion of lemming winter nests with signs of reproduction also varied considerably among years and habitats (Supplementary Material S1 Fig. S6). We found evidence that proportion of nests with reproduction was 1.2 and 3.3 times higher in riparian compared to mesic and wet habitats, respectively, (Table 1b) and was twice higher in collared than in brown lemmings (β = 0.69, CI = [0.45, 0.94]; Fig. 3a). Results from the years when data from the shrub habitat was available revealed the same trends and suggested a low reproductive activity in shrubs although no significant difference was found, probably due to small sample size (Table 1b, Fig. 3b). However, reproductive rate was still 2 times higher in collared lemmings than in brown lemmings in this time series (β = 0.70, CI = [0.11, 1.29]).

When examining the influence of snow parameters on lemming reproduction, we found support for an increase in reproduction with greater snow depth and a reduction in reproduction as maximal height of depth hoar increased (Table 2b, Fig. 4 a-b, Supplementary Material S1 Table S2). A relationship with the inverse of snow depth was preferred over a linear one (ΔAICc = 1.8), suggesting a positive effect on reproduction until ~60 cm, after which there was little change (Fig. 4a). We also found evidence of a decrease in reproduction with an increase in density and thermal conductivity of the basal snow layer (Fig. 4c-d; Table 2b).

Variation in snow properties across habitats

The density of the basal snow layer was highest in riparian habitat, lowest in shrub, and intermediate in the other two habitats (Fig. 5a; coefficients of all models are shown in Supplementary Material S1 Table S3). We did not find any evidence of differences in the snow thermal conductivity (keff) across habitats, but we note a high residual variability of the random variable (SDresidual = 8.77), suggesting large variations within years (Fig. 5b). Snow SSA was 1.2 to 1.4 times lower in the riparian habitat compared with the others (Fig. 5c). Snow depth differed between all habitats and, as expected, it was 2.8 to 3.1 times deeper in riparian habitat compared to the other three habitats (Fig. 5d). Maximal height of depth hoar was lowest in the wet habitat and highest in riparian and shrub habitats (Fig. 5e). Finally, snow onset dates in the riparian habitat were almost two weeks earlier than in the wet habitat (Fig. 5f).

We found strong evidence that temperature of the basal snow layer was warmest in the riparian habitat, coldest in wet and intermediate in shrub CI = [1.7,8.32 S3, Fig. S8). When correcting for the slight difference in air temperature between the shrub and other habitats (1 o C warmer in shrubs on average), we found similar results (Supplementary Material S1 Fig. S9).

The principal component analysis (PCA) revealed an association between deep snow cover, thick depth hoar and high density of the basal layer, for which riparian habitat globally had higher scores than the other three habitats (Supplementary Material S2, Table S4, Fig. S11). There was also an association between thermal conductivity and SSA, and scores of the riparian habitat on this axis indicated a negative association with these variables (see Supplementary Material S2 for more details on the results).

Discussion

Our results show that habitat use by lemmings in winter is not random and that snow properties can explain some of the observed patterns. Reproductive activity also differs among habitats and can also be partly explained by spatial variations in snow properties. Overall, riparian habitat was the one with the highest quality in winter as revealed by the highest nest density and reproductive rate, followed by shrub being intermediate, and wetland and mesic habitats being the lowest quality. These patterns were broadly similar for both lemming species despite some differences, and they did not change with population density.

On the importance of a warm and soft snowpack

The first weeks of the winter are thought to be a critical period for lemmings since temperature drops rapidly and the thin, heterogenous snowpack does not provide a significant protection from cold temperatures. Thus, lemming may seek habitats where snow accumulates first, as is the case in riparian and shrub habitats. The riparian habitat also had the deepest snowpack and the thickest depth hoar, two snow characteristics known to offer good insulating properties [START_REF] Zhang | Influence of the seasonal snow cover on the ground thermal regime: An overview[END_REF]. Indeed, we observed more stable and warmer temperature in the basal snow layer of riparian compared with other habitats (up to 17 °C warmer compared to wetland in the coldest months). These conditions should decrease thermoregulatory costs for lemmings [START_REF] Chappell | Thermal energetics and thermoregulatory costs of small arctic mammals[END_REF]) and could therefore explain why riparian habitat had the highest use and reproductive activity. These observations are consistent with previous studies highlighting the importance of snow depth in lemming habitat use [START_REF] Batzli | Habitat use by lemming near Barrow , Alaska[END_REF], Duchesne et al. 2011a[START_REF] Schmidt | Long -term patterns in winter habitat selection, breeding and predation in a densityfluctuating, high Arctic lemming population[END_REF][START_REF] Beckerath | Long-term monitoring reveals topographical features and vegetation explain winter habitat use of an arctic rodent[END_REF]) but also in other mammal species like wolverine (Gulo gulo; [START_REF] Glass | Spatiotemporally variable snow properties drive habitat use of an Arctic mesopredator[END_REF] or brown bear (Ursus arctos; [START_REF] Sorum | Den-site characteristics and selection by brown bears (Ursus arctos) in the central Brooks Range of Alaska[END_REF].

Use of the riparian habitat could nonetheless come with a cost as we found evidence for higher snow density in its basal layer compared to other habitats. Since thermal conductivity and density are correlated [START_REF] Domine | Linking the effective thermal conductivity of snow to its shear strength and density[END_REF], we were surprised that thermal conductivity was not higher in the riparian habitat, but this might be related to technical difficulties when taking this measurement manually in friable depth hoar layers [START_REF] Fourteau | On the use of heated needle probes for measuring snow thermal conductivity[END_REF]. We did find smaller SSA values in this habitat, suggesting a greater extent of indurated depth hoar formed in meltfreeze layers, which is consistent with our observations. As snow accumulates early in the riparian habitat, this makes it more susceptible to melting episodes in early winter when temperatures are fluctuating around 0 °C. Furthermore, the vertical temperature gradient is reduced in deeper snowpacks, which slows down depth hoar formation and leads to a denser and less developed depth hoar layer [START_REF] Marbouty | An experimental study of temperature-gradient metamorphism[END_REF]. Therefore, the dense and relatively hard basal layer found in deep snow could explain the reduced reproduction that we observed.

For lemmings, staying inside their nests in winter offers an energetic advantage [START_REF] Chappell | Thermal energetics and thermoregulatory costs of small arctic mammals[END_REF], so their foraging trips should be as quick and efficient as possible. Digging in hard and dense snow increases their energy expenditure but also decreases their digging speed [START_REF] Poirier | Snow hardness impacts intranivean locomotion of arctic small mammals[END_REF]), which may increase the time spent outside the nest and thus thermoregulatory costs. Our results suggest that the complex spatial variability in snow properties may impose tradeoffs on lemmings because even if deep snow areas provide the most favorable thermal environment, they may not offer the best conditions to minimize travel costs. Ultimately, this may reduce their energy available for reproduction. In the boreal forest, a similar tradeoff was observed in the Pacific marten (Martes caurina), which prefers habitats with deep snow as a shelter against the cold despite the increased locomotion costs associated with this type of habitat [START_REF] Martin | Forest structure and snow depth alter the movement patterns and subsequent expenditures of a forest carnivore, the Pacific marten[END_REF].

Despite the known preference of brown lemmings for wet habitats in summer [START_REF] Batzli | Habitat use by lemming near Barrow , Alaska[END_REF], our study shows that it is the least used habitat in winter, likely due to its thin snowpack and high daily temperature fluctuations (Duchesne et al. 2011a). Nonetheless, this summer preference could explain the trend for a greater use of this habitat by brown compared to collared lemmings in winter. In addition, considering that lemmings favor horizontal movements over vertical ones when digging tunnels in the snow [START_REF] Poirier | What guides lemmings movements through the snowpack?[END_REF], the presence of sharp mounds surrounding polygons in the wet habitat (Fig. S1) may impede lemming movements.

Our results suggest that the shrub is a relatively good habitat for lemmings in winter. Shrubs not only favor early and relatively deep snow accumulation, but they also reduce densification due to overburden compaction, therefore providing optimal conditions for the development of soft depth hoar where movement is facilitated [START_REF] Domine | The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime[END_REF]). In addition, food resources such as willow buds may encourage lemmings to use this habitat [START_REF] Soininen | Highly overlapping winter diet in two sympatric lemming species revealed by DNA metabarcoding[END_REF][START_REF] Fauteux | Feeding preference of brown lemmings (Lemmus trimucronatus) for plant parts of Arctic willow (Salix arctica)[END_REF]). We did not find evidence for a greater reproductive activity in this habitat, but this may be related to the relatively small number of nests found during the 3 years of sampling the shrub habitat. Nonetheless, the favorable snow conditions offered by this habitat will strongly depend on the weather at the beginning of the winter. Indeed, if early winter warm spells lead to the formation of a hard melt-freeze snow layer, this will preclude snow from drifting, preventing the accumulation of a deep snow cover around shrubs [START_REF] Barrere | Snowmelt events in autumn can reduce or cancel the soil warming effect of snow-vegetation interactions in the arctic[END_REF]. It is worth noting that the favorable subnivean temperature found in this habitat was enhanced by the slightly warmer air temperature in that area of the study site compared to the other sites due to the reduced cool katabatic flow [START_REF] Domine | Permafrost cooled in winter by thermal bridging through snow-covered shrub branches[END_REF]. Moreover, because this habitat is spatially limited in the High Arctic, it is difficult to extrapolate our results in this habitat to other systems.

It is also worth noting that our study focused only on snow properties but other factors such as vegetation or soil characteristics could also influence lemming distribution and reproduction in winter (Duchesne et al. 2011a). Nonetheless, experimental studies showed that snow per se is a strong factor affecting lemming distribution independently of other landscape features [START_REF] Reid | Lemming winter habitat choice : a snow-fencing experiment[END_REF]).

Interspecific differences and density dependence

We did not find any difference in habitat use between the two lemming species during winter despite known differences in their summer habitat use [START_REF] Batzli | Habitat use by lemming near Barrow , Alaska[END_REF]. One reason could be that both lemmings have a relatively similar winter diet at our study site dominated by dicotyledons, which are mainly found in riparian, shrub and mesic habitats [START_REF] Soininen | Highly overlapping winter diet in two sympatric lemming species revealed by DNA metabarcoding[END_REF]. We found a higher winter reproductive activity in collared lemmings compared to brown lemmings, independently of the habitat. Collared lemmings are thought to be better adapted to the extreme conditions of the High Arctic [START_REF] Fuller | High-arctic lemmings, Dicrostonyx groenlandicus. II Demography[END_REF] and to be more efficient at digging in hard snow compared to brown lemmings [START_REF] Poirier | Snow hardness impacts intranivean locomotion of arctic small mammals[END_REF], which might help them save energy for winter reproduction. Considering that winter reproduction is generally a key determinant of lemming outbreaks [START_REF] Millar | On reproduction in lemmings[END_REF], Fauteux et la. 2015), it is surprising that it is the brown and not the collared lemming that reaches the highest densities at our study site in those years [START_REF] Gauthier | Long-term monitoring at multiple trophic levels suggests heterogeneity in responses to climate change in the Canadian Arctic tundra[END_REF]. This suggests that other demographic factors such as mortality may limit population growth in collared lemmings. This species is known to be more vulnerable to avian predators compared to brown lemmings during the summer [START_REF] Seyer | Resource partitioning among avian predators of the Arctic tundra[END_REF], and a high reproductive activity may increase predation rate in winter nests [START_REF] Fauteux | Density-dependent demography and movements in a cyclic brown lemming population[END_REF].

Contrary to [START_REF] Schmidt | Long -term patterns in winter habitat selection, breeding and predation in a densityfluctuating, high Arctic lemming population[END_REF], we did not find evidence of density dependence in habitat use. This suggests that when conditions are favorable for the development of a high-quality basal snow layer (Domine et al. 2018), there may be enough refuges and resources to support both lemming species even in peak years at our study site. Nonetheless, brown lemmings are thought to be dominant over collared lemmings during competitive interactions and can force them to use sub-optimal habitats [START_REF] Morris | Measuring the ghost of competition: Insights from density-dependent habitat selection on the co-existence and dynamics of lemmings[END_REF]. Brown lemmings may sometimes take over collared lemming nests and even perform infanticide if young are present (M. Poirier, pers. obs.).

Therefore, considering the large overlap in habitat use in winter, we cannot exclude that more subtle differences due to competitive interactions may still exist between the two species. We should also bear in mind that this study focuses on habitats that act as shelter for lemmings (i.e., where they build their nests), which could differ in some cases from the habitats where they forage.

Conclusion

Understanding the winter ecology of arctic species is a huge challenge, especially those living within the snowpack. Our study differs from previous ones by directly assessing the links between snow properties and lemming habitat use instead of relying on proxies of snow conditions, like topographic features or habitat classes [START_REF] Vaillant | Spatial distribution in Norwegian lemming Lemmus lemmus in relation to the phase of the cycle[END_REF][START_REF] Schmidt | Long -term patterns in winter habitat selection, breeding and predation in a densityfluctuating, high Arctic lemming population[END_REF][START_REF] Beckerath | Long-term monitoring reveals topographical features and vegetation explain winter habitat use of an arctic rodent[END_REF]. We showed that spatial variation in snow properties can affect lemming habitat use and reproduction, and thus could affect their demography, as previously shown for annual variations in snow conditions [START_REF] Bilodeau | The effect of snow cover on lemming population cycles in the Canadian High Arctic[END_REF], Domine et al. 2018b). Moreover, our results highlight the complex links between snow properties and the ecology of animals living in this medium as we found that habitats offering the most favorable microclimate (deep snow) may not offer the best conditions for reproduction due to the presence of harder snow. The Arctic snowpack is expected to be strongly affected by climate change as rain-onsnow and melt-freeze events should increase at higher latitudes and cause a hardening of the basal snow layer, especially in critical habitats where snow accumulates early in winter [START_REF] Hansen | Warmer and wetter winters: Characteristics and implications of an extreme weather event in the High Arctic[END_REF]. Arctic regions exposed to oceanic currents and with milder and wetter climate are especially exposed to such events [START_REF] Boonstra | Why do the boreal forest ecosystems of Northwestern Europe differ from those of Western North America?[END_REF]). In the worst cases, basal ice crust could encapsulate vegetation and deprive herbivores from their food [START_REF] Berteaux | Effects of changing permafrost and snow conditions on tundra wildlife: critical places and times[END_REF], leading to massive mortality as was observed in reindeer populations [START_REF] Stien | Congruent responses to weather variability in high arctic herbivores[END_REF][START_REF] Langlois | Detection of rain-on-snow (ROS) events and ice layer formation using passive microwave radiometry: A context for Peary caribou habitat in the Canadian Arctic[END_REF][START_REF] Dolant | Assessment of the Barren Ground caribou die-off during winter 2015-2016 using passive microwave observations[END_REF]. On the other hand, warming of the Arctic will lead to shrub expansion [START_REF] Tape | The evidence for shrub expansion in Northern Alaska and the Pan-Arctic[END_REF], which could provide favorable habitats for small mammals in winter.

Researchers should aim to monitor snow properties in the long-term but also at a fine spatial scale to better assess the determinants of the population dynamics of Arctic mammals.

Figure legends

  ]; βshrub-wet = 4.53, CI = [1.21, 7.85]; n = 72; Fig 6a, Supplementary Material S1 Table S3, Fig S7). Similarly, we found that daily temperature fluctuations were lowest in riparian habitat, highest in wet and intermediate in shrub, suggesting a more insulative snowpack in riparian and shrub (βriparian-shrub = -1.70, CI = [-2.03, -1.37]; βshrub-wet = -1.04, CI = [-1.36, -0.71]; n = 144; Fig 6b, Supplementary Material S1 Table

Fig. 1 Fig. 3 Fig. 4

 134 Fig. 1 Lemming winter nest density in three different habitats (2007 -2019) for (a) brown

Fig. 5

 5 Fig. 5 Snow properties of four winter habitats used by lemmings at Bylot Island, 2014 -2019.

Fig. 1

 1 Fig. 1

Fig. 2

 2 Fig. 2

Fig. 3

 3 Fig. 3

  Fig. 4

Fig. 6

 6 Fig. 6

Acknowledgement:

We thank Gabriel Bergeron, James Akpaleeapik, Madelaine Proulx, Qalaapik Enookolo, Jaimie Vincent and all the other people involved in the long-term monitoring of lemming winter nests in the field. We also thank Mathieu Barrère and Marianne Valcourt for their help with snow sampling. We thank Marie-Christine Cadieux for her precious help regarding database management. We thank Christophe Kinnard for sharing Reconyx data with us and Denis Sarrazin for his technical support with automated weather stations. We also thank David Bolduc for his help with the map in appendix. Funding: This work was funded by Sentinel North program of the Canada First Research Excellence Fund, the Fonds de recherche du Québec -Nature et technologies, the Natural Sciences and Engineering Research Council of Canada, the Polar Continental Shelf Program of Natural Resources Canada, the Fondation de l'Université Laval and the W. Garfield Weston Foundation.

contribution from GG, FD, DF. GG and FD cosupervised the project and obtained funding.

Availability of data:

The lemming data analysed in this study are available from the NordicanaD repository (https://doi.org/10.5885/45400AW-9891BD76704C4CE2) and the other datasets are available from the corresponding author.

Conflict of Interest:

The authors declare that they have no conflict of interest.

Table 1 (a) Coefficients of the models examining the influence of three (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) or four habitats (2017)(2018)(2019) on lemming nest density (ln-transformed). (b) Coefficients of the models examining the influence of three (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)