Study area
We studied elk from four Colorado elk herds, the Avalanche Creek (Colorado Parks and Wildlife Data Analysis Unit [DAU] E-15), Bear’s Ears (DAU E-2), Trinchera (DAU E-33), and Uncompahgre Plateau (DAU E-20) herds that all wintered in Colorado, USA (Fig. 1). In 2021, population estimates for each herd were: Avalanche Creek ~ 4800, Bear’s Ears ~ 18,400, Trinchera ~ 14,000, and Uncompahgre Plateau ~ 12,500. Elk from the Bear’s Ears herd infrequently crossed into Wyoming, while elk from the Trinchera herd often crossed into New Mexico. Similar to other partially migratory elk populations, migrants in these systems tend to overwinter at low elevations and venture to higher elevation ranges for the summer. Resident elk, by contrast, remain year-round on one range or regularly commute between two ranges.
Elevation in the study area ranges from ~ 1500 m at lower elevations to > 4000 m in mountainous regions. Vegetation communities are diverse and largely determined by altitude and aspect32. Low elevations (below 2000 m) harbor pinyon pine (Pinus edulis)-Utah juniper (Juniperus osteosperma) forests, plains grasslands, sagebrush steppe, high desert, desert basins, and river valleys. Moderate elevations (2000–3500 m) host montane shrubs, oak (Quercus gambelii), sagebrush (Artemisia spp.), aspen (Populus tremuloides) or mixed aspen/conifer, Engelmann spruce (Picea engelmannii), and fir (Abies spp.) communities, while high altitudes (> 3500 m) are characterized by alpine systems or talus slopes. The study area is a mosaic of public and private land ownership.
Animal capture and monitoring
We captured yearling and adult female elk on winter ranges by net-gunning and darting them from a helicopter during January-March 2017–2021. We radio-collared elk with Global Positioning System (GPS) collars programmed to acquire a location every 1, 2, or 4 h (G2110E2, G5-2D, Advanced Telemetry Systems, Isanti, MN, USA) and monitored most individuals for 13–18 months. We collared 104 elk from the Avalanche Creek herd (2019–2021; ~ 2% of the estimated population), 82 elk from the Bear’s Ears herd (2019–2021; < 1% of the estimated population), 113 elk from the Trinchera herd (2017–2021; ~ 1% of the estimated population) and 158 elk from the Uncompahgre Plateau herd (2017–2021; ~ 1% of the estimated population). Additionally, we outfitted most pregnant females with vaginal implant transmitters (VITs; M3930U, Advanced Telemetry Systems, Isanti, MN, USA) during capture. See Appendix Table A1 for annual collaring totals. In all of our analyses, each elk was only included once, i.e., each elk only contributed one elk-year of data to each analysis. This study was carried out in accordance with ARRIVE guidelines33. All capture, handling, and collaring of elk were completed in accordance with applicable guidelines and regulations and were approved by the Colorado Parks and Wildlife Animal Care and Use Committee (protocol IDs: 02–2017, 01–2019, 02–2019, 01–2020, 03–2020).
Identifying migrants and residents
We followed the methods of Bunnefeld et al. (2011) and used time-series of net-squared displacement (NSD) methods and positional elevation to classify movement strategies of individual elk34. We categorized elk as either migrants, residents, commuters, or dispersers. We classified individuals that exhibited discontinuous season-specific space use with movement events between winter range and summer range as migrants. We classified individuals that inhabited one range year-round or that exhibited regular year-round movement between proximate ranges (generally along a short-distance elevational gradient) as residents. We classified individuals that demonstrated long-distance movement away from their initial location without return as dispersers. Individuals falling outside of all established categories were considered ‘ambiguous’ (Appendix Figure B1). We determined seasonal space use of each elk using the timing of their migrations, producing estimates of departure and arrival dates by fitting movement models with the R package MigrateR35 (see Appendix A for more details).
Variable extraction
For covariate extraction purposes, we included elk-years of migratory individuals with clear spring departure and arrival dates, complete spring migrations in their data, and at least two weeks of space use on both sides of spring migration events (i.e. pre-departure period in the winter range and post-arrival period in the summer range). We selected a two-week period to characterize the spatiotemporal conditions of forage immediately prior to departing from winter range and immediately after arriving on summer range. Elk-years with incomplete tracking data were included in migration metrics calculations (Table B1) when possible but excluded from other statistical analyses. Because collaring in the Uncompahgre Plateau herd occurred after spring migration had already begun, we could not estimate departure timing in the first year of monitoring, and instead used data from the second year. We extracted spatiotemporal covariates over three periods: 1) the two weeks prior to departure from winter range, 2) the spring migration, and 3) the two weeks following arrival at summer range (hereafter pre-departure, during migration, and post-arrival respectively). We used occurrence distributions calculated from continuous-time movement models to delineate space use in each period36. We determined the distance between winter and summer ranges by calculating the distance (in kilometers) between the centroids of the two-week periods of space use on either side of the spring migration.
We calculated instantaneous rate of green-up (IRG) metrics from the MODIS09Q1 Version 6 MODIS/Terra Reflectance product (250 m, 8-day resolution) by fitting a double sigmoid curve to 8-day Normalized Difference Vegetation Index15. From these curves, we extracted both interpolated daily IRG values and year-specific estimates of the ordinal day of peak annual IRG per pixel11,15. We calculated the mean per-location IRG optimality and the mean ordinal day of peak annual IRG over the spatial area encountered during the pre-departure, migration, and post-arrival periods. We calculated the instantaneous rate of green-up optimality for a given pixel and elk location as:
$${OPT}_{y,t,p}=1-(\text{max}\left({IRG}_{y,p}\right)-{IRG}_{y,t,p})$$
where OPT is optimality, y is the year, t is the ordinal day of the year, p is the pixel containing a given relocation coordinate, and IRG is the instantaneous rate of green-up value. The index theoretically ranges from 0 to 1, with 1 indicating perfect optimality, i.e., no difference between contemporaneous forage quality and the maximum forage quality value at that location for the year. This metric better reflects the current condition of forage in comparison to its maximum potential than the commonly used Days-From-Peak metric, because locations can have identical peak days, but green-up at different rates18. Additionally, we calculated the standard deviation in peak IRG date to represent variability in date of peak green-up along the migration route as a way of characterizing the greenscape18. We also extracted snow covariates from the same MODIS product. We calculated the per-pixel snow-free DOY representing the day-of-year where the absence of snow was first detected in an 8-day resolution time series for a given pixel3.
We estimated calving dates for elk in our analysis when data were available (n = 128). We estimated calving dates primarily using the expulsion date of vaginal implant transmitters (VITs; M3930U, Advanced Telemetry Systems, Isanti, MN, USA) inserted into pregnant females during capture. Occasionally, VITs failed or were expelled prematurely. In these cases, we attempted to use GPS collar location data from pregnant females to identify parturition events. If we verified a parturition event by capturing a calf or finding a birth site at a GPS cluster, we used the timestamp of the pregnant female’s first GPS location at that cluster as the estimated parturition date. We did not have calving date estimates for all elk-years used in our analyses because we were not able to identify a calving date for every elk-year, and because we sometimes used GPS datasets from elk-years in which elk no longer had VITs. As a result, analyses involving calving date sample sizes were reduced (Avalanche Creek n = 38, Bear’s Ears n = 59, Trinchera n = 12). Uncompahgre Plateau elk were excluded because we could not relate calving date to spring migration timing in this herd.
Departure and duration models
Our two response variables were departure date from the winter range and duration of the spring migration. Departure date was defined as the day-of-year (DOY) from January 1st of the focal year, and duration was defined as the number of days between departure from winter range and arrival at summer range. We estimated departure date if net-squared displacement data was asymptotic prior to displacements indicative of a migratory event (representing seasonal home-ranging behavior) to avoid including elk collared during their migration.
To generate models investigating potential determinants of spring departure date and the duration of migration, we used eight covariates as independent variables (see section above about how they were extracted or calculated). Productivity covariates included the mean peak IRG date pre-departure, the mean peak IRG date during migration, the mean peak IRG date post-arrival, and the standard deviation of peak IRG date along the migration route. Snow metrics included the mean snow-free date pre-departure, the mean snow-free date during migration, and the mean snow-free date post-arrival. Additionally, we included the distance between seasonal ranges as a covariate.
We used multiple linear regression to fit herd-specific models. We did not use a generalized linear mixed model framework that would have pooled data across herds due to convergence issues when including multiple random slopes to estimate herd-specific responses37. Our candidate models included every linear combination of each covariate with no more than four covariates per model. Models were limited to four covariates to avoid overfitting given our herd sample sizes. In the Trinchera herd we lowered this to a maximum of two covariates per model due to a lower sample size of migratory individuals relative to other herds (n = 17 elk). In each model we excluded covariates with Pearson’s correlation coefficient ≥ 0.7 (i.e. correlated variables were not included in the same model). We selected variables within correlated pairs by comparing univariate model fits and excluding the covariate whose model resulted in significantly lower explanatory power as determined by Akaike Information Criterion for small sample sizes (AICc)38. Where univariate models were equally parsimonious (ΔAICc < 2), we ran each version of the uncorrelated candidate model (i.e. models retaining one or the other correlated variable) and included them both in model selection. We fit each candidate model and selected the top model using AICc. A null, intercept-only, model was also included in the list of candidate models of elk migration phenology. We chose as our top model the most parsimonious model where ΔAICc was within two of the top-ranking model to identify the most important potential drivers (Appendix Table B2). We removed covariates with a Variance Inflation Factor ≥ 5 and removed influential outliers from the data39. We identified potential outliers using diagnostic plots and removed outliers if they were both abnormal and influential; n = 10 elk-years removed from models (Bears Ears n = 6, Uncompahgre Plateau n = 4). For models incorporating calving date as a covariate, our fitting and selection methods were identical to those used for our primary models but excluded the Uncompahgre Plateau herd, given that we could not relate calving date estimates to spring migratory timing in this herd (Appendix C).
Variation in optimality of access to forage
To explore inter-individual variation in exposure to IRG pre-, during-, and post-migration and identify potential sub-strategies, we used Gaussian-mixture model clustering from the R package mclust40. We used the mean estimated optimalities per individual during each respective period (pre-, during, and post-migration) as input covariates. Top models and the optimal number of clusters were automatically selected via Bayesian information criterion (BIC)38. We performed a chi-squared goodness of fit test using herd and optimality cluster as categorical variables to investigate the null hypothesis that there was no association between herd and optimality sub-strategy, i.e., we compared the observed frequency of herd in each cluster to a null expectation that each herd was divided evenly between clusters.
