Denning phenology and reproductive success of wolves in response to climate signals

We compiled GPS location data from 388 individuals in eight wolf populations across western Canada and Alaska (figure 1) (Longitude: −154.6 to −112.7°, Latitude: 51.5 to 67.8°; table 2). Wolves were captured following standard animal care protocols defined by affiliate university or government agencies and released with GPS collars programmed to acquire fixes over a range of intervals from 15 min to 24 h. Each population was monitored from 1 to 12 years between 2000 and 2017.

Zoom In Zoom Out Reset image size

The movements and site-fidelity patterns of resident pack members can signal pup-rearing activities from March through August (Alfredéen 2006, Tsunoda et al 2009). These patterns emerge from food provisioning and other social interactions that regularly draw pack members back to pup-rearing sites that serve as social centers throughout the summer (i.e. dens and rendezvous sites) (Ballard et al 1991, Alfredéen 2006, Ciucci and Mech 2006). Thus, the movements of non-reproductive individuals can be as informative as reproductive individuals (Tsunoda et al 2009).

We classified movements indicative of pup rearing using a three-step process. For the first step, we identified sites with a high-intensity of use by clustering collared animal locations in space (⩽100 m) and time (⩽7 d) (R package rASF (Mahoney and Young 2017)). This process is performed iteratively over all locations within an animal's truncated time series, producing convex hulls for each set of points that meet the criteria. We then merged all overlapping convex hulls (in space and time) into a single cluster and removed any clusters with fewer than eight locations (or five for 24 h fix interval datasets). For the second step, we smoothed the same movement time series using a median filter and an overlapping, 4 d moving window to dampen the effect of large movements. We flagged 4 d periods with median daily displacements less than 200 m and overlaid the output on our cluster data. For the final step, we visually inspected movement time series to evaluate whether cluster fidelity persisted after frequent offsite forays (i.e. provisioning of offspring) and to identify any influential gaps in location data that could affect estimates for den initiation dates (figure 2). We estimated parturition as the initial date from the earliest denning cluster observed during the reported period for wolf parturition (March–June; Mech and Boitani 2010). Occasionally, movement data from a reproductive female expressed gaps in a location time series that lasted approximately 4 d to 8 d, an indication of GPS satellite signal occlusion while in underground dens (Joly et al 2018). Visual inspection of the time series helped to identify these gaps so that den initiation dates (i.e. parturition dates) could be adjusted to when the signal was first lost. However, in cases where denning initiation were unclear, particularly for known non-breeding individuals, the movement data were removed from further consideration.

Zoom In Zoom Out Reset image size

We define reproductive success as packs with one or more pups at the end of August. To estimate success, we evaluated the movement time series for each individual after a denning event had been identified. If pup-rearing movements were evident through the end of August using the same methods above (e.g. high-fidelity clusters and low daily median displacement through August), we flagged these packs as reproductively successful. Although these activities often extended well into September or October, we chose August as a more conservative date given the variability in when pups begin to consistently travel with natal packs. We validated success by comparing our estimates to visual observations of pups with packs during the autumn or winter for a portion of our dataset for which these observations were made (Denali National Park and Preserve, Banff National Park, Jasper National Park, and West Athabasca River study areas).

In addition, because we performed these assessments for each individual independent of pack membership, we used pack members as a form of validation in our estimates for both den initiation and reproductive success. However, for the analyses below, we used only one estimate of reproductive timing and success per pack per year, prioritizing dates derived from reproductive females followed by individuals with the highest quality data (i.e. highest fix retention or sampling rate).

We aggregated seasonal weather metrics using Daymet (v3; Thornton et al 2018), a meteorological product that contains daily estimates for minimum and maximum temperatures, total precipitation, and snow water equivalent (SWE) on a global 1 km grid. To characterize vegetation dynamics during the most photosynthetically active periods (i.e. growing season), we estimated Normalized Difference Vegetation Index using 8 d surface reflectance derived from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS; MOD09Q1 data product, ORNL DAAC 2017). We post-processed NDVI data using the program TIMESAT (Jönsson and Eklundh 2004), which masked cloud-covered pixels, smoothed NDVI time series, and estimated phenological metrics such as start of growing season (SOS), length of growing season (LOS), and time-integrated NDVI (tiNDVI). We defined growing season start and end dates by when pixels passed the 10% and 90% thresholds for mean NDVI amplitude as measured across a seasonal time series. We estimated snow disappearance date (SDD), or the first snow free day after a minimum of 3 d with snow cover (assessed backwards in time to capture the end of the snow-covered season), using the MODIS normalized snow difference index (MOD10A1; Hall et al 2006) and the Google Earth Engine API. We also included PDO (http://jisao.washington.edu/pdo; Mantua and Hare 2002) and Arctic Oscillation indices (AO) (NOAA National Weather Service Climate Prediction Center; www.cpc.ncep.noaa.gov) during January of the reproductive year and annual means from the previous year (effectively, a lag of one year; figure 3(b)). See table 3 for a complete description of each covariate.

Zoom In Zoom Out Reset image size

To define local domains for weather conditions, we calculated seasonal wolf home ranges using 95% isopleths from fixed-kernel density estimates (Sheather and Jones 1991). We defined three biological seasons: summer (pup-rearing: April–August), autumn (ungulate rut: September–November), and winter (wolf pair formation and breeding: January–March). We used seasonal home ranges to extract median climate statistics, and home range centroids to estimate latitude for each individual. In the absence of wolf movement data during any non-denning seasons, summer home ranges summaries were used instead given the strong seasonal site fidelity exhibited by wolves (Mech and Boitani 2010).

Prior to the analyses, we assessed collinearity across all covariates and removed one or more covariates that contributed to correlation coefficients greater than 0.7. In cases where correlation occurred, we retained the metric that best represented the underlying, hypothesized biological mechanisms. We also evaluated temporal trends in median SOS and denning date during the years for which we had wolf denning data by fitting linear mixed effects model with either SOS or denning date as the response variable, a single continuous effect of year, and random intercepts for study system. We estimated shifts in dates based on predictions derived from models with significant effects of year (95% confidence interval for year did not overlap zero).

Next, we evaluated inter-seasonal correlations (figure 1). We defined population domains by their minimum convex polygons (MCP; Calenge 2006) using all collared individuals within a population. For each population domain, we estimated median SOS, seasonal temperatures, cumulative precipitation, and SWE across all years within the study (2000–2017). We then generated population-specific Pearson's correlation coefficients by pairing SOS with climate metrics from the previous winter, autumn, and summer during each year, thus producing annual estimates of correlation for each climate metric by population across all years.

We used a time-to-event model, Cox proportional hazard regression (R package survival, Therneau 2015), to examine the effects of climate factors on denning phenology. Denning phenology was measured as the number of calendar days since January 1st and the initiation date of each documented denning event. The baseline 'hazard' in this context reflected the daily probability of denning from January through the middle of June. We derived robust, Huber sandwich variance estimates to account for non-independence in the timing of denning for packs with more than one denning event across multiple years (i.e. specified 'cluster(PackID)' in the model formulation; Therneau 2015). We evaluated the influence of seasonal (previous summer, autumn, and winter) minimum and maximum temperatures, previous summer and autumn cumulative precipitation, mean daily SWE, SDD, SOS, LOS, tiNDVI, latitude, annual PDO, and annual AO on the timing of wolf den initiation.

We assessed climatic factors affecting reproductive success using mixed logistic regression with success (1) and failure (0) as a binary response (R package lme4; Bates 2010). We included random intercepts for pack nested within population and evaluated the influence of synchrony with spring onset (i.e. the difference between denning date and start of growing season), annual PDO (previous year), seasonal weather conditions, and time-integrated tiNDVI. In contrast to the phenology analyses, summer conditions were based on the current reproductive year rather than from the previous summer. In addition to a linear effect of spring onset, we considered a non-linear second-order polynomial and a linear piecewise response to synchrony with spring onset (i.e. using a single knot at 0 representing perfect synchrony) to account for possible differences in success whether litters were born before or after spring onset. We also included summer home range area (km²) as an indirect measure of home range quality, because smaller home ranges are often associated with higher resource productivity (Duncan et al 2015).

We standardized all covariates by centering on the mean and dividing by one standard deviation prior to modeling (Gelman et al 2014). We ran all-possible additive models using the R package MuMIn (Barton 2009). Models were ranked by Akaike Information Criterion corrected for small sample size (AICc; Burnham and Anderson 2002). We generated model-averaged coefficient estimates and confidence intervals based on unconditional standard errors following Burnham and Anderson (2002). Model-averaged coefficient estimates were derived across all models, but the estimate for any particular covariate was conditional on its presence within a model (i.e. coefficients were not fixed at 0 when absent). We determined variable importance using parameter weights (Burnham and Anderson 2002) and confidence interval overlap with zero (85% level; Arnold 2010). We also evaluated both model sets for uninformative parameters using personally authored R code (sensu Leroux 2019), which we retained when deriving parameter weights (to achieve covariate balance across model sets; Burnham and Anderson 2002) and model-averaged estimates in an all possible models context (Arnold 2010).

We assessed model goodness-of-fit by using the Cox proportional hazard concordance statistic for the denning phenology model set (Therneau 2015) and conditional pseudo-R² for the reproductive success model set (Nakagawa et al 2017). Concordance is a measure of agreement between observed and predicted values commonly used in time-to-event models (Therneau 2015). This is accomplished by ranking the 'risk' scores (exponentiated linear predictor in a hazard model) of a pack at each observed denning event relative to all packs who have yet to den. Values of zero correspond to the lowest 'risk' of denning relative to the remaining packs (i.e. poor predictive ability); whereas, values of one correspond to the highest 'risk' of denning (i.e. perfect predictive ability). Concordance is then estimated as the weighted average of these rankings across all denning events, with values approaching one indicative of stronger predictive accuracy.

We compiled movement data from 388 wolves and identified 227 possible dens associated with 106 packs across western Canada and Alaska between 2000 and 2017. Of these, we classified reproductive success for 186 reproductive events in those packs with sufficient movement data (i.e. data from denning through the end of August and ⩾1 fix per day). Although rare, we treated multiple litters as a single event based on the timing of the first den and any pup-rearing movements through August as a single measure of success. Validations comparing our movement-based predictions to aerial and ground observations of denning indicated that we successfully identified 100% of the known denning events (n = 146 known dens). We had two or more collared individuals for 41 of these dens, permitting us to calculate variation in estimated den date across individuals of different sex and breeding status. The median difference was 1 d (µ = 1.95, SE = 0.31, range = 0 d to 8 d).

Of 153 confirmed aerial or ground observations of recruitment (i.e. packs with or without pups after August 31st of each year), 143 cases matched our predictions of success (93.5%). Of the ten misidentified, six were estimated as failures but observed to be successful and three were estimated as successful but observed to be failures. Thus, error in recruitment classification was relatively minor without, we believe, inducing any systematic bias associated with the 'type' of error due to near equal representation. We corrected the ten misidentified reproductive success estimates to reduce the overall classification error rate in our data set.

Autumn and winter mean temperatures were negatively correlated with start of growing season (Autumn: r_mean = −0.50, SE = 0.07; Winter: r_mean = −0.40, SE = 0.05) and did not show a trend over the 18-year period for which we had data (2000–2017) (figure 4). Summer temperature was also negatively correlated with SOS (r_mean = −0.57, SE = 0.09) but demonstrated a slight trend toward weaker correlation in recent years. Summer and winter precipitation, as well as winter SWE, showed slight but highly variable correlations with SOS and no apparent temporal trend (figure 4). However, autumn precipitation exhibited a temporal trend, switching from positive to negative correlation during the same 18 year period.

Zoom In Zoom Out Reset image size

The estimated median denning date was May 4 (SD = 13.6 d), on average 14.7 d prior to the start of season. SOS advanced an average of 14.2 d from 2000 to 2017, whereas the median denning date did not change (figures 3(a) and S1 (available online at stacks.iop.org/ERL/15/125001/mmedia)). Denning was initiated earlier at lower latitudes, but variability in the timing of denning occurred among all populations. Latitude was strongly correlated with seasonal temperatures (e.g. Winter T_max: r = −0.77), as was SOS and SDD (r = 0.83). Because effects of latitude and SDD were largely accounted for by the other climate variables (figure S2), latitude and SDD were removed from further consideration in this analysis. Further, we limit our discussion to seasonal maximum temperatures because of consistent model interpretations and generally improved model fit over minimum temperatures. Approximately 18.1% of phenology models expressed some degree of uninformative covariates (Leroux 2019). When present, these covariates consisted of LOS (5.5%), SOS (2.3%), autumn precipitation (4.3%), summer precipitation (3.7%), and autumn temperature (2.3%), which were consistent with relatively low parameter weights for these covariates (table 5). Models with evidence of uninformative covariates did not appear in our top models (table 4).

The best model (ΔAICc relative to null = −43.52; Concordance = 0.67, SE = 0.02) for denning phenology included winter and summer maximum temperatures, tiNDVI, annual PDO, and autumn precipitation in order of support via parameter weights (tables 4 and 5). In general, denning occurred earlier in regions or years with warmer temperatures during summer and winter seasons preceding parturition. However, denning was delayed following years with higher tiNDVI and PDO, as well as with higher precipitation in autumn and winter (table 5).Importantly, our results may indicate that wolf reproductive phenology was responding to photoperiod, or latitude as a proxy, rather than winter temperatures (Asa and Valdespino 1998). Although statistically intractable due to a strong correlation between latitude and winter temperature, we included an additive effect of latitude within our best reproductive phenology model to evaluate how the two metrics were confounded. Notably, the model maintained significant effects for winter (β = 0.25; SE = 0.13) and summer temperatures (β = 0.27; SE = 0.09), but estimated a non-significant effect for latitude (β = −0.20; SE = 0.15). Though qualitative, this result would indicate there was some redundancy between temperature and latitude (i.e. photoperiod), but that seasonal temperatures were capturing additional variation independent of latitude; therefore, temperatures were more informative drivers of phenology at both regional and local population levels.

The best model for reproductive success included autumn precipitation, summer maximum temperatures during pup-rearing, annual PDO, and home range area (in order of support via parameter weights) (ΔAICc = −15.85 relative to the null model; $R_{}^2 = 0.22$; table 5). Model-averaged coefficient estimates indicated pup recruitment through the end of August increased significantly with positive PDO during the previous year and declined with increased autumn precipitation, increased summer maximum temperatures during the current pup-rearing season, and increased home range area (tables 4 and 5; figure 5). We found no support for a linear effect, second-order polynomial, or a linear piecewise response to synchrony with spring onset (with a knot at perfect synchrony). The piecewise analysis allowed testing for differences in reproductive success whether animals denned early or late relative to the mean difference between den date and SOS. The confidence intervals for both piecewise coefficients overlapped 0 (β_<0 = −0.23, CI₉₅ = −1.05 to 0.59; β_>0 = 0.21, CI₉₅ = −0.42 to 0.84), indicating no difference in the influence of phenological match on reproductive success between cases where denning occurred early or late relative to SOS. Approximately 35.1% of success models expressed some degree of uninformative covariates (Leroux 2019). When present, these covariates consisted of annual AO (13.0%), autumn temperature (10.2%), winter temperature (5.4%), match with SOS (3.2%), summer precipitation (2.6%), and tiNDVI (0.7%), which were also consistent with relatively low parameter weights for these covariates (table 5). Models with uninformative covariates did not appear in our top models (table 4).

Zoom In Zoom Out Reset image size