Statistical Analyses:
All statistical analyses were conducted in R v.4.0.3 (R Core Team, 2020) using packages climwin v1.2.3 (Bailey & van de Pol, 2016; van de Polet al., 2016), lme4 v1.1.26 (Bates et al., 2015), asreml-R v4.1 (Butler, 2021), ggplot2 v.3.3.3 (Wickham, 2016), tidyverse v.1.3.1 (Wickham et al., 2019), nadiv v.2.17.1 (Wolak, 2012), lubridate v.1.7.10 (Grolemund & Wickham, 2011), gridExtra v.2.3 (Auguie, 2017), ggeffects v.1.1.2 (Lüdecke, 2018), lmerTest v3.1.3 (Kuznetsova et al., 2017), and car 3.0.10 (Fox & Weisberg, 2019).
One of the key assumptions in a marmot’s life history is that they reproduce immediately following their emergence. We wanted to investigate this assumption by estimating the relationship between female emergence from hibernation and pup emergence date from weaning in a given year. To do so, we used a restricted dataset from the years 2003-2017 where we had all female emergence dates from hibernation and pup emergence from weaning in each year. We then ran a linear mixed model with the package lme4 (Bates et al., 2015) using pup emergence date as our response variable. We used female emergence date as the only fixed effect. Year and female identity were fitted as random effects. Since we were expecting a 1:1 ratio between pup and female emergence date, we also used a t-test to compare the slope of the female emergence date to 1.
To determine during which phenological window temperature and snowpack had the strongest association with pup emergence date, we used the R package climwin (Bailey & van de Pol, 2016; van de Pol et al.,2016). For both environmental variables, we fit a model in climwin with pup emergence date as our response variable and either average temperature or snowpack as our independent variable. Our specified reference day was June 1st, and the starting date of the window varied from June 1st to November 13th (200 days before June 1st) the previous year. We allowed any length of window from 1 to 200 days.
To investigate whether climate change was impacting the timing of reproduction in the yellow-bellied marmot, we fitted a univariate animal model of pup emergence date using the asreml-R package (Butler, 2021). For this model, we restricted the dataset to include those females of known age, with a known mass in June, who reproduced that year. Fixed effects included mother’s age, valley (up- or down-valley), average spring snowpack, average spring temperature, and the mother’s mass in June. Random effects were the year, permanent environment (see Kruuk & Hadfield, 2007), additive genetic, and colony effects. Year was added to control for inter-annual variation in conditions experienced. Permanent environment effect was added to control for any inter-individual variation in pup emergence date not due to genetic effects. Additive genetic effects were added to estimate the amount of variation in the phenotype associated with additive genetic variation. Colony was added to control for potential micro-environmental differences between them. Significance of the random effects was determined using a log-likelihood ratio test. Starting from our full model, random effects were dropped one at a time and the log-likelihood ratio of each model were compared.
To investigate inter-individual variation in the plasticity of the timing of reproduction, we extended the previous animal model by adding random slopes for average spring snowpack and temperature in two separate models using asreml-R (Butler, 2021). We tested the significance of the random slope terms and the random intercept by comparing the log-likelihood of the models with and without each of these effects. We did not include additive genetic effects in this model.
We estimated the existence of selection (directional and/or stabilizing) on pup emergence date. To do this, we ran generalized linear mixed models with three different fitness proxies in lme4 (Bates et al., 2015). For each litter, we used the total number of pups surviving to one year old, the proportion of pups surviving to one year old, and the total number of pups. We used the same model structure for the three models. Fixed effects were the linear and quadratic orthogonal polynomials of pup emergence date, mother’s age, mother’s mass in June, average spring snowpack, average spring temperature, and valley. Random effects were female identity, colony, and year. We used a Poisson distribution for number of pups surviving to one year old and the number of pups in the litter. We used a Binomial distribution for the proportion of pups surviving to one year old.
For all models, all continuous variables fit as fixed effects were scaled to have a variance of one and a mean of zero.

Results: