Density-Dependent Foraging Behaviors on Sympatric Winter Ranges in a Partially Migratory Elk Population

Introduction

Density dependence plays a central role in understanding a variety of ecological processes including population dynamics and community organization. Density influences habitat selection and movements by animals to determine species distributions at both small and large scales (Owen-Smith et al., 2010; Almeida et al., 2015). Partial migration, where individuals follow resident or migratory tactics of movement (Dingle and Drake, 2007), is hypothesized to result from density-dependent trade-offs between costs and benefits influenced by phenotype, individual state, or the behavior of conspecifics in the population (Berg et al., 2019). Migratory and resident tactics therefore may be maintained by differential density-dependent regulation of vital rates that must demographically balance each other over the long term (Kaitala et al., 1993; Hebblewhite and Merrill, 2011). Further, there is limited evidence that migratory tactics are genetically fixed, and recent data show that at least in some populations ungulates switch between migratory tactics (Berg et al., 2019), and that switching may be density dependent (Eggeman et al., 2016). How density-dependent changes in behaviors of large herbivores might contribute to switching and more broadly the maintenance of partial migration has not been well-studied.

In temperate montane systems, migratory ungulates typically move to areas of high-quality resources in summer (Mysterud et al., 2011; Middleton et al., 2013; Monteith et al., 2018), where they meet reproductive demands and gain body reserves to improve survival and reproduction during the resource-limited winter (Cook et al., 2004, 2013; Monteith et al., 2014). In such systems, most studies focus on exposure to forage resources and predation risk when ungulates are on separate ranges to explain density-dependent maintenance of partial migration, and assume individuals derive equal costs and benefits when together on their sympatric range (e.g., Kaitala et al., 1993; Ball et al., 2001; Hebblewhite and Merrill, 2009). However, foraging behaviors of migrant and resident ungulates may differ when on sympatric ranges, which may offset fitness advantages of migration (Robinson and Merrill, 2013; Found and St. Clair, 2016). Such differences may result from differential exposure to external conditions or social grouping, and the demographic consequence may be exacerbated by conspecific density, which in most montane systems typically is higher in winter than in summer because of snow constraints (Cagnacci et al., 2011). Further, if migration is flexible with individual variation attributed to year-to-year condition-dependence in a facultative manner (Nelson, 1995; Fieberg et al., 2008; Grovenburg et al., 2011; Monteith et al., 2011, but see Sawyer et al., 2019), then a better understanding of the forage-predation trade-offs migrants and residents make as population densities change may contribute to our understanding of the maintenance of partial migration.

Under density-dependent habitat selection, non-territorial animals are expected to distribute themselves in a manner that achieves equal fitness within habitats, which is altered by the density of conspecifics present (Fretwell and Lucas, 1970; Rosenzweig, 1981; Morris, 2003). Studies have shown that ungulates change their habitat selection across densities (McLoughlin et al., 2006; Perez-Barberia et al., 2013; Mansson et al., 2017), but habitat-selection analyses alone may not reveal the behavioral mechanisms behind trade-offs in the selection process (Bastille-Rousseau et al., 2010). Strong selection of resource-rich patches is expected to emerge because foragers exhibit area-restricted search whereby they decrease the speed and increase the tortuosity of movements that results in high patch-residence times, as compared with moving directly through resource-poor patches (Benhamou, 1992; Anderson et al., 2008; Avgar et al., 2011). If resource depletion increases with conspecific abundance, herbivores may be forced to explore the same areas more intensively to obtain the necessary resources, which would result in higher concentration and predictability of individuals in high-quality areas (Almeida et al., 2015). Alternatively, individuals may move more to explore larger areas, leaving a patch after a shorter time (Jiang and Hudson, 1993; Shipley and Spalinger, 1995; Searle et al., 2005), making them less predictable but more likely to encounter predators (Daly et al., 1990; McKenzie et al., 2012). Individuals are expected to return when patch resources recover or depletion in other patches increases the value of the patch (Barraquand and Benhamou, 2008; Van Moorter et al., 2009; Seidel and Boyce, 2015). Concomitantly, in gregarious species, high rates of interaction and conflicts of direction in movement of individuals may limit group speeds (Fortin et al., 2009; Pays et al., 2012; Sigaud et al., 2017), which in turn may put constraints on exploring the environment, reducing foraging intake, or result in group fission (Lardy et al., 2016). Herbivores make compromises between foraging and safety (Lima and Dill, 1990; Verdolin, 2006; Visscher et al., 2017), leading to the generality that risky foraging decreases (Abrams, 1993; Brown, 1999). To escape predation, prey alter movement strategies to reduce predator encounters in time and space (Kie and Bowyer, 1999). Van Moorter et al. (2016) reported the strength of selection of resources by female moose (Alces alces) was associated with high variation in both residence and return times reflecting resource heterogeneity, but their study took place under low risk from natural predators. In contrast, residence times of North American elk (Cervus canadensis) in Wisconsin were unrelated to selection but reflected frequent, unpredictable returns to preferred areas on the landscape, which was hypothesized to be a strategy that reduced predation (Anderson et al., 2008).

Group size and cohesion reflect how foragers balance limiting factors (Visscher et al., 2017). Indeed, by associating with large groups, individuals can dilute the per capita probability of mortality (Hamilton, 1971), provide “many eyes” for predator detection (Dehn, 1990; Brown and Kotlar, 2004), or reduce foraging costs from vigilance, at least under encounter-limited foraging (Fortin et al., 2004; Robinson and Merrill, 2013). But foraging in large groups can have other consequences, such as increasing conspecific interference or rates of depletion resulting in increased movement rates (Molvar and Bowyer, 1994; Kausrud et al., 2006; Mobaek et al., 2012). These trade-offs may depend on the social relationships that exist among animals within a group (Weckerly, 1999; Millspaugh et al., 2004; Vander Wal et al., 2013, 2014). Therefore, changes in abundance that alter grouping patterns are likely to alter the outcomes of trade-offs between reducing predation risk and accessing forage (Fortin et al., 2009; Pays et al., 2012).

In this paper, we evaluated the assumption that migrant and resident elk make similar trade-offs between foraging opportunities and safety from wolves on the sympatric winter range of the Ya Ha Tinda adjacent to Banff National Park (BNP) (Figure 1) in response to ∼70% population decline over a 14-year period (Hebblewhite et al., 2018). Previously, we have reported migrant elk are exposed to higher forage quality on summer ranges consistent with the forage maturation hypothesis (Hebblewhite and Merrill, 2008), and indicating migrants may enter winter in better condition. This assumption is supported by migrants having higher fecal nitrogen, pregnancy rates, and calf weights (Hebblewhite and Merrill, 2011). At the same time, on winter ranges we observed differences in foraging behaviors of migrant and resident that constrained synchronizing food processing with vigilance, which we hypothesized could alter foraging benefits between migratory tactics particularly under encounter-limited foraging (Robinson and Merrill, 2013). If changes in winter density differentially influences forage-predation trade-offs of migrant and resident elk, differences in winter foraging could contribute to demographic balancing among migratory tactics (Hebblewhite and Merrill, 2011; Eggeman et al., 2016).

To test hypotheses concerning density-dependent foraging, we focused on diurnal trade-offs in selection for abundant forage areas under predation risk, and on changes in social cohesion and group sizes because elk are most active during the day in winter (Supplementary Figure S2). We assumed diurnal periods corresponded to time spent foraging. We predicted elk would exhibit density-dependent habitat selection, but we expected the extent of response to be less for residents because they needed to compensate for low forage quality in summer (Hebblewhite and Merrill, 2009). Because residents maintained larger, more cohesive groups than migrants (Robinson et al., 2010; Hebblewhite et al., 2018), we also expected the strength of selection for high-forage areas by residents would be associated with longer residence times and high return times. In contrast, because migrant elk maintained less coherent groups that were less familiar with risk from predators including humans (Robinson and Merrill, 2013), we predicted greater trade-offs between high-forage areas and safety and increased movements in risky areas. Elucidating behavioral mechanisms associated with changes in density on populations of partially migratory ungulates in winter is key to broadening our understanding of how forage trade-offs on seasonal ranges interact and contribute to condition-dependent switching and the maintenance of partial migration (Berg et al., 2019).

Materials and Methods

Study Area

Elk winter range is located on the eastern slopes of the Canadian Rocky Mountains adjacent to BNP, and overlaps with Parks Canada’s 4,000-ha Ya Ha Tinda horse ranch (Hebblewhite et al., 2006 for details). The study area is montane temperate, with temperatures in the area averaging 9°C during May–September, and −4.1°C during winter. Precipitation averages 319 mm during summer (May–September) with an average of 157 cm of snowfall. During the study, resident elk remained on the grassland year-round, whereas migrants moved either west into Banff National Park or east to multiple-use forests (Eggeman et al., 2016). During this study, the partially migratory elk population declined from ∼1,400 animals in 2002–2003 to ∼450 in 2015–2016. The leading cause of the decline was high predation by gray wolves (Canis lupus) on adult females (Hebblewhite et al., 2006, 2018), combined with high predation by grizzly (Ursus arctos) and black bears (Ursus americanus) on neonatal elk (Berg, 2019). Switching among migratory tactics occurs but is relatively low (15%/yr; Eggeman et al., 2016). Elk were the most abundant ungulate during our study, but other ungulates included moose (Alces alces), mule deer (Odocoileus hemionus), white-tailed deer (O. virginianus), and bighorn sheep (Ovis canadensis). The main predators of elk in winter when our observations occurred were wolves (5.7–8.9 wolves/1,000 km²; Hebblewhite, 2006) and First Nations peoples in winter. Cougars (Puma concolor) also were active during winter and thought to be increasing in this region (Knopff et al., 2014).

Elk Population Size

Elk population size was surveyed every winter as trend counts within the same area when migrant and resident elk were on their sympatric winter range (Eggeman et al., 2016). Aerial surveys were conducted by Alberta Fish and Wildlife (2002–2010) or Parks Canada (2011–2016) using rotary-wing or fixed-wing aircraft every winter (January–February after heavy snowfall to maximize sightability) except for the winter of 2012–2013 and 2015–2016 when we used maximum counts of elk on winter grounds (Supplementary Table S1). Aerial counts were adjusted for sightability by a 15% increase based on previous research (Hebblewhite, 2000). In most years, these numbers were repeatedly verified from the ground.

Elk Capture and Monitoring

Adult female elk were captured and fitted with Global Positioning System (GPS) collars (Lotek GPS 3300, 4400, 7000; LOTEK Inc., Aurora, ON, Canada; n = 92) or very high frequency (VHF) collars (LMRT-4, LOTEK Inc., Aurora, ON, Canada; n = 147) between February 2002 (winter 2001–2002) and March 2015 (winter 2014–2015) on the winter range using either corral trapping, helicopter net gunning or darting from horseback (Animal Care Protocols: University of Alberta 353212, 611812, 000624, University of Montana AUP 004-16, Parks Canada BAN-2014-16756). We collected GPS fixes at a variety of intervals, but rarefied data to a consistent 2-h schedule. Location error (∼34 m) and fix-rate bias (<10%) were low enough to avoid bias in the selection models (Hebblewhite and Merrill, 2007; Hebblewhite and Merrill, 2008). Although all collared elk migrated back to the Ya Ha Tinda to winter, not all elk arrived by 1 November (Supplementary Figure S1). We included only those locations of migrant elk that arrived on Ya Ha Tinda by 1 November and only female elk that had at least 100 locations during the winter (155 elk winters from 92 adult female elk). We omitted potentially erroneous GPS fixes (<1%, including unreasonably long-steps) and GPS fixes that were not preceded by two consecutive GPS fixes were removed in the selection and movement analyses because three consecutive fixes were required to calculate turn angles. We classified individuals as migrant or resident each year using all locations from the previous spring and summer based on the net squared displacement (NSD) method (Bunnefeld et al., 2011; Borger and Fryxell, 2012) in the R package MigrateR (Spitz et al., 2017), combined with post-hoc spatial rules and visual confirmation by plotting movements of individual animals (Eggeman et al., 2016 for details). Data are available on Movebank (Eggeman et al., 2016).

We defined winter as 1 November to 31 March because 95% of migrants completed their fall migrations and were on the winter range by 1 November, and 100% of the elk were still on the winter range on 31 March (Eggeman et al., 2016). The diurnal period was defined as 0600–1900 h when elk most actively foraged based on their movements (Supplementary Figure S2). Selection and movement analyses were based on 155 elk winters of data from 92 unique GPS-collared individuals (Supplementary Table S2). The average number of days tracked for each elk-winter was 58.4 (SD = 42.9; range 9–151) days with a mean of 608.0 (SD = 450.5; range 107–1,806) locations.

Movement-Independent Resource Selection

We used integrated Step Selection Analysis (iSSA) to estimate selection for areas varying in forage and predation risk for individual elk while controlling for movement rates based on fixed-time step lengths and directional persistence based on turn angles. The iSSA method is an extension of the step selection function approach (Fortin et al., 2005; Thurfjell et al., 2014) that accounts for animal movement; failure to adequately control for the movement process may lead to biased resource-selection estimates (Forester et al., 2009). The iSSA simultaneously estimates movement and resource-selection parameters allowing for effects of environmental variables on selection processes to be distinguished from movement (Avgar et al., 2016). We used the R package amt (Avgar et al., 2016) to conduct the iSSA using a binomial link. We fit a gamma distribution to observed step lengths and a Von Mises distribution to observed turn angles using Maximum Likelihood methods. We randomly drew 10 available steps per observed step using these distributions. Covariate values for forage biomass and wolf predation risk were then extracted for the endpoints of all observed and available steps.

Wolf predation risk (hereafter, predation risk) was quantified using a diurnal, resource-selection function (RSF) for winter developed from GPS and VHF telemetry data from five wolf packs (n = 30 wolves) in the Ya Ha Tinda from 2002 to 2004 and weighted by the kill rate/pack/day (see Hebblewhite, 2006 for details). The wolf RSF was modeled as a function of landcover, elevation (m), slope (%), aspect class (north, south, flat), distance to forest edges (km), and distance to human activity (km). We extended the predation risk to 2005–2016 by incorporating changes over time in annual wolf pack size, land cover, and distance to edge but used the same RSF model to predict selection (Berg, 2019). We determined that predation risk using this approach was correlated to wolf-scat surveys conducted later in the study in summer (i.e., 2013–2016; Spilker, 2018), and assumed similar correspondence in winter.

Herbaceous (graminoids and forbs) biomass (g/m²) was predicted across Ya Ha Tinda during 2002–2004 from empirical GLM models based on landcover type, MODIS greenness, and elevation that varied annually (Hebblewhite, 2006; Hebblewhite and Merrill, 2008). For 2005–2015, we used the predicted forage biomass layer from 2002 to 2004 that corresponded most closely in terms of summer precipitation, but updated estimates of forage biomass in burns (prescribed and wildfires: Banff Fire Database 2015, Parks Canada) and clearcuts (Sundre Forest Products, 2015) using a forage biomass-successional model (Berg, 2019 for details). We assumed late summer spatial patterns in forage biomass corresponded to forage spatial patterns in winter (Hebblewhite, 2006).

We fit conditional logistic regression models for each elk in each winter (individual elk-winters) using the function amt:fit_issf in the R package amt and summarized selection by migratory tactic. We included movement variables (ln steplength and cosine turn angle), forage biomass, predation risk, and their interaction in all models because we were explicitly testing for differences in the same variables between migrants and residents rather than finding the best model for each individual. Both forage biomass and predation risk variables were standardized by subtracting the mean and dividing by the standard deviation to allow for comparisons of effect size among individuals. To test our hypotheses that elk selection changed with elk population size, we modeled the iSSA-derived coefficients of individual elk for forage biomass, predation risk, and their interaction separately using migratory tactic (migrant as reference) and annual elk abundance (counts) of elk on the Ya Ha Tinda winter range. We used Generalized Linear Mixed-effect Models (GLMMs) with elk identity as a random intercept assuming a Gaussian distribution, and weighted the model coefficients for each elk-winter by 1/variance of the coefficient (Gillies et al., 2006).

Daily Movement, Resident and Return Times

We used elk diurnal locations (0600–1900; Supplementary Figure S2) to calculate daily step lengths (m), residence times (hrs), and return times (hrs) for all GPS-collared elk used in the resource selection analyses. We based step lengths on the straight-line distance between 2-h relocations of GPS-collared elk. Residence time was the duration of time spent within the vicinity (800-m radius) of any diurnal location and was estimated using the package adehabitatLT (Barraquand and Benhamou, 2008). We used a circle of 800-m radius as the spatial threshold because that distance approximated the median distance (m) moved during the active period, and outputs based on radii from 400 to 1,600 m were correlated (r > 0.30, p < 0.005). We calculated return times as the amount of time (hrs) it took to revisit an area within a 200-m radius of a previous GPS location after 48 h using the R package recurse (Bracis et al., 2018). We determined these threshold values for return times produced qualitatively similar relationships with elk population size (Supplementary Figure S3). The mean residence and return times were calculated for all the locations within the movement path of each individual during winter, except elk that never returned to a previously used area; these elk were not included in the return-time analysis. This procedure reduced the number of elk used in the analysis by one migrant and three resident-elk winters.

To associate step length, residence time, and return time to environmental conditions, we related each metric to forage biomass or wolf predation risk at the start-point location, migration status (migrant = reference), and their interaction using a GLMM with elk identity and year as random effects assuming a Gaussian distribution. We used Akaike Information Criterion (AIC) to compare the fit of linear and non-linear (quadratic) fit of the models. To assess our hypotheses that overall movement characteristics differentially changed with elk population size between migratory tactics, we used GLMMs, assuming a Gaussian distribution, to model each movement metric including elk population size, migratory tactic, and their interaction as covariates, elk identity as a random intercept.

Elk Associations and Group Sizes

Because data for elk group size were not associated with the GPS locations of elk, group size could not be explicitly integrated into selection and movement analyses (sensu Berger et al., 2015). Instead, we determined group sizes visually from the ground at an interval of approximately 7–10 days during most winters using VHF-telemetry to locate both VHF- and GPS-collared elk (Hebblewhite and Merrill, 2011). We defined a group associated with a collared individual as one to n individuals that interacted with each other, headed in the same direction, or used a confined area at the time of observation. We determined elk group size by scanning with binoculars or telescope and using a clicker-counter to count individuals. We assessed the change in the distribution of group sizes with elk population size with quartile regression (Proffitt et al., 2012; Brennan et al., 2015), and differences in the median of group sizes between periods of high (>600 elk) and low (<600) elk population size within migratory status using a median test. We used a break point of 600 because the population had declined by approximately one-half its size at the beginning of the study (Supplementary Table S1).

We assessed association among classes of dyads of GPS-collared migrant and resident elk using the coefficient of sociality (C_s, Kenward et al., 1993) within the R package wildlifeDI (Long et al., 2014). C_s varies between −1 and 1 and is calculated for dyads of GPS-collared elk as:

where d_O is the mean Euclidean distance between simultaneous locations of an elk dyad, and d_E is the mean Euclidean distance of n² permutations of the simultaneous fixes of a dyad (Kenward et al., 1993). For this analysis, we considered locations that were within 1 h as “simultaneous” because collars took fixes every 2-h and some collars did not fix on the same hour. To assess our hypothesis that migrant elk were more associated with resident elk at low elk abundance, we used GLMs to relate C_S to elk population size and dyad type (resident-migrant, migrant-migrant, and resident-migrant; migrant-migrant as reference) and their interaction.

Results

Movement-Independent Resource Selection

Selection for high-forage areas was consistent across nearly all (99%) elk, whereas most (76%) elk selected for rather than against areas of high predation risk (Table 1), likely because forage biomass and predation risk were positively correlated (r = 0.44, p < 0.001, n = 1,000). Nevertheless, 75% of the elk showed a trade-off in forage to avoid wolf predation by reducing their selection for high-forage areas under high risk (i.e., negative F × W interaction; Table 1). There were no main or interactive effects of migratory tactic in how elk selected for areas (p > 0.52).

As elk population size increased, the strength of elk selection for areas of high forage biomass decreased (β = −0.00020, SE = 0.000040, p < 0.001; Figure 2A). In contrast, elk increased their selection for risky areas (β = 0.00016, SE = 0.000032, p < 0.001; Figure 2B) as elk population size increased, but they also exhibited a weaker trade-off (Figure 2C) for high-forage areas as predation risk increased (β = −0.00057, SE = 0.000027, p < 0.036; Supplementary Figure S4). Migrant and resident elk showed similar changes in selection with elk population size (p > 0.34).

Diurnal Movement Rates, Residence, and Return Times

Consistent with strong selection for forage biomass, elk reduced step lengths (β = −44.66, SE = 5.61, p = 0.001; Figure 3A), increased residence time (β = 0.94, SE = 0.24, p = 0.001; Figure 3C), and had shorter return times (β = −38.05, SE = 3.57, p = 0.001; Figure 3E) to areas of high forage biomass during the foraging period. Residents altered step lengths less, residence time less, and return time more in response to forage biomass than migrants (Figure 3A, Supplementary Table S3, F × M interaction). In contrast, elk exhibited more complicated and nonlinear responses to increasing predation risk (Supplementary Table S3). At low predation risk, elk had long step lengths (Figure 3B), low residence time (Figure 3D), and long return times (Figure 3F) that initially declined (step length and return times) or increased (residence time) as predation risk increased, which is similar to responses to increasing forage biomass. Nonetheless, at very high predation risk, the slopes of the relationships changed or tended to level off (Figure 3 and Supplementary Table S3). Overall, migrant and residents responded similarly in step lengths with predation risk, whereas residence time and return times were higher for residents than migrant (Figure 3 and Supplementary Table S3).

Mean step lengths of elk increased (β = 0.121, SE = 0.042, p = 0.004; Figure 4A) and residence time of elk decreased (β = −0.11, SE = 0.001, p = 0.001; Figure 4B) with elk population size, with no interaction with migratory tactic (p > 0.26). As elk population size changed, there was no consistent change in return times for elk following either movement tactic (Figure 4C and Supplementary Table S4).

Elk Group Sizes and Associations

We counted the number of individuals associated with each marked elk from 1 November to 30 March in 59 to 717 groups/yr across the 14-year period. Only the upper 75th quartile of group sizes decreased with decreasing elk population size (β = 0.134, SE = 0.005, p < 0.001), with group sizes of resident elk decreasing more than migrants (β = 32.0, SE = 4.15, p < 0.001). Median group size with which marked elk were associated was smaller at high (>600 elk) than at low elk population size (<600) in both residents (201 vs 275, median test: X² = 49.9, p < 0.001) and migrants (171 vs 253, X² = 51.91, p < 0.001). A difference in median group size existed between migrant and resident elk only when elk population size was low (253, X² = 7.35, p = 0.007; Figure 5).

Estimates of the C_S were derived for 14 ± 22 (mean ± SD) resident-resident, 14 ± 22 migrant-migrant, and 38 ± 64 migrant-resident GPS-collared elk dyads/year across the 14-year period. Overall, migrant-migrant and migrant-resident dyads had similar coefficients of sociality (C_S, p = 0.26), which were mostly higher than resident-resident dyads (p = 0.04, Supplementary Table S5). As density decreased, coefficients of sociality of migrant-migrant and migrant-resident dyads increased more than for resident-resident dyads (Figure 6).

Discussion

Contrary to our predictions, we documented only minor differences in the selection and movement responses of migrant and resident elk to forage biomass and predation risk on sympatric winter ranges. Both migrant and resident elk consistently showed selection for areas of high forage biomass during the daytime when elk spend most of their foraging time (Ensing et al., 2014; Kohl et al., 2018; Supplementary Figure S2). The strong selection for areas of high forage biomass emerged from longer residence times and shorter return times (Almeida et al., 2015). Shorter return times to areas of high forage biomass compared with low biomass areas have been reported previously during the growing season in other elk populations (Wolf et al., 2009; Seidel and Boyce, 2015), but not during winter. In contrast, elk showed a weaker and less consistent response to spatially avoiding predation risk. Indeed, we noted positive selection for areas of high predation risk overall, which we attribute to the correlation between forage biomass and predation risk in our and other ecosystems (e.g., Yellowstone; Kohl et al., 2018). Only at the very highest levels of predation risk did we observe elk moving more quickly through risky areas and not returning to them as soon, which was reflected in reduced selection for areas of high forage biomass as predation risk increased. The same trade-off in resource selection was reported for resident elk in summer (Hebblewhite and Merrill, 2009).

We expected stronger selection responses by elk to predation risk because wolves are the leading cause of mortality for adult females in this system (Hebblewhite et al., 2018). Several factors may contribute to the moderate response. First, the extent of response may reflect the frequency of reactive vs proactive responses to the risk of predation (Creel et al., 2014; Moll et al., 2017; Cusack et al., 2020). Because we used a RSF-based predation risk, that variable reflected primarily habitat conditions associated with wolves (i.e., risky places), where elk are expected to proactively make trade-offs in selection or alter their movements because of perceived vulnerability (Kauffman et al., 2007; Creel, 2018; Kohl et al., 2018). Based on the strong spatial pattern of vigilance of elk at Ya Ha Tinda (Robinson and Merrill, 2013), we conclude that elk at the Ya Ha Tinda perceived predation risk, but they avoided only the riskiest places during the day by moving through them quickly. Cusack et al. (2020) reported that elk in Yellowstone National Park did not exhibit proactive avoidance of predation risk, which they attributed to wolves moving frequently throughout the winter range of elk such that elk had imprecise knowledge of their locations. In contrast, Creel et al. (2005) showed that where elk are directly confronted by wolves, elk shifted out of preferred grazing sites in open meadows into the cover of coniferous woodland. We did not have simultaneous locations of wolves and elk to address fine-scale temporal responses, but Robinson and Merrill (2013) also reported elk vigilance was not as closely related to wolf presence or sign in the area as to where elk were located, specifically relative to distance to timber and human infrastructure. We also did not have comparable data to assess predation risk from cougars, but cougars were detected only four times as compared with wolves 98 times on 44 remote cameras on the winter range in 2017–2018 (Flowers, 2019).

If elk were frequently subject to direct threats of predation by the wolves during diurnal foraging at Ya Ha Tinda, we might have observed stronger trade-offs for safe areas. Fewer encounters may result from reduced movement of wolves during the day, but our analysis included crepuscular periods, when wolf activity has been reported to be highest (Krawchuk, 2014). Nevertheless, wolves generally remained on the periphery of the open grasslands during the daytime in winter, which we have attributed to wolves’ avoidance of human activity (Hebblewhite and Merrill, 2009; Robinson et al., 2010). These so-called “human-shield” areas that seem capable of allowing herbivores to exploit resources in safe areas are widespread across ungulate systems (Ogutu et al., 2005; Berger, 2007; Rogala et al., 2011; Whittington et al., 2019). Thus, elk exhibited less reactively to risk, but avoided “risky places” by moving through them quickly. In either instance, Ya Ha Tinda elk could more fully exploit forage resources across the productive grassland during the daytime. In contrast, by shifting their distribution to use areas near human infrastructure at night when wolves move into the grassland (Hebblewhite and Merrill, 2009; Robinson et al., 2010), night-time responses are more consistent with the “risky times” hypothesis; similar responses of elk have been documented in Yellowstone to temporal variation in both wolf and cougar predation risk (Kohl et al., 2019).

Second, elk at Ya Ha Tinda maintained large aggregations (median > 200 elk across all years), which allowed them to avoid predation even if directly confronted by wolves. Creel et al. (2008) argued that aggregation is a type of proactive response to predation. Large prey aggregations are common in open habitats (Clutton-Brock et al., 1987; Pays et al., 2007; Gower et al., 2009), which are expected to reduce the per capita risk to an individual as well as the number of overall encounters (Huggard, 1993; Hebblewhite and Pletscher, 2002; McLellan et al., 2010). Because homogeneously distributed resources are thought to facilitate decisions in movement directions and speed, Sueur et al. (2011) hypothesized that groups should be more cohesive in such landscapes. We argue that the maintenance of large group sizes and the availability of a human-mediated refuge from wolves has allowed an increasing number of resident elk to remain on the Ya Ha Tinda year-round while the population has declined (Hebblewhite and Merrill, 2011).

At the same time, we documented few differences in selection and movement responses of migrant and resident elk as density changed. Our prediction that resident elk would exhibit weaker density-dependent habitat selection to compensate for lower forage quality in summer was not supported. Instead, migrant and resident elk similarly diminished their selection for areas of high forage biomass as elk abundance increased (Morris, 2003; Fortin et al., 2008). Thus, our results add to the growing empirical evidence for density-dependent selection in large herbivores (Clutton-Brock et al., 1987; McLoughlin et al., 2006; Perez-Barberia et al., 2013; van Beest et al., 2014). The inverse relationship between residence times and elk population size indicated elk did not intensify search efforts within the same area as population size increased, but rather that interference or exploitive competition forced them to leave patches sooner and explore other areas to meet their requirements (Almeida et al., 2015). Fortin et al. (2009) argued that exploitative competition among foraging bison was insufficient to explain reduced time spent in grassland patches of even large groups of bison because vegetation removed during a foraging bout was <6%. We suggest short residence times during winters at high elk population size reflected cumulative forage depletion because over-winter forage use was as high as 94% during the period of highest elk abundance (McInenly, 2003). Patterns of snowfall may complicate this interpretation, because low snowfall years coincided mostly with years of high elk population size. With reduction in snow, when access to forage is greater, we would expect longer residence time (rather than travel), but this is not the case. We suggest that short residence time is because of forage depletion.

If forage depletion were the major driver of increased foraging movements at high elk population size, we also would have expected the weakest forage-predation risk trade-offs because individuals would become increasingly food stressed (Peacor, 2003). Instead, the trade-off of high-forage areas for safety intensified as elk population size increased (Figure 1C and Supplementary Figure S4). One explanation is that individual elk were associated with a broader range of group sizes at high elk population size (Figure 5). At high population size, large groups may be associated with increasing probability of group fission (Body et al., 2015) because of intra-group competition and reduced coordination in travel (Pépin and Gerard, 2008; Fortin et al., 2009; Pays et al., 2012). This may account for relatively more small groups at high elk population size. We cannot rule out a relatively higher detection bias in counting small groups vs large groups in the field, but because groups were tracked by finding collared animals, we still would have located small groups even if we missed some individuals within the group. Instead, shifts in the portions of large and small group sizes with elk population size likely reflects the nonlinear relationship between density and group sizes described for gregarious ungulates across a number of systems (Mansson et al., 2017).

As the population of elk declined, we were surprized that the median group size increased rather than decreased (McLellan et al., 2010). Because wolf predation on elk is not density dependent in this system and wolf abundance did not decline (Hebblewhite et al., 2018), per capita predation risk from wolves may have increased as elk abundance declined leading to stronger gregariousness. Maintaining relatively large groups under high predation risk reflects an optimal group size in terms of trade-offs in group foraging constraints, reducing predator encounters, and maximizing dilution effects (Focardi and Pecchioli, 2005). Associated with the increase in group size, we observed an increase in spatial cohesion among elk, particularly between migrants and between migrants and residents. In contrast, cohesion remained more constant in resident animals across populations sizes. This outcome might be expected because residents maintain large groups year-round whereas migrant groups are smaller and more dispersed in summer (Hebblewhite et al., 2006; Hebblewhite and Merrill, 2011). More empirical studies linking sociality, habitat selection, and population density are needed to understand the generality of what we observed in this study (Webber and Vander Wal, 2018).

Implications for Maintaining Partial Migration

In the face of widespread decline of migratory behavior in ungulates, there is growing emphasis in understanding the behavioral mechanisms that maintain partial migration to facilitate conservation efforts (Nicholson et al., 1997; White et al., 2007; Rickbeil et al., 2019; Sawyer et al., 2019). Partial migration is thought to be maintained by density-dependent fitness balancing between strategies at the population-level or by switching at the individual level (Lundberg, 1988; Kaitala et al., 1993). Previously, we concluded that maintenance of partial migration of elk at the Ya Ha Tinda was a result of demographic balancing of vital rates despite very different resource selection strategies in summer to escape predation and find forage. At the time, it was unknown whether year-round residents would have the foraging advantages once on their sympatric winter ranges. Subsequently, we have garnered little support for differences in space use or foraging behaviors by elk in winter (Robinson et al., 2010, this study) that might offset the nutritional carry-over effects that migrant elk accrue in summer. Thus, we expected that as density declined, residents should exhibit the most improvement in reproduction, which did not occur (Hebblewhite et al., 2018). Instead, we observed density-independent predation by wolves, which we have argued has kept elk well below ecological carrying capacity (K) such that little improvement in recruitment would be expected as density has declined (Hebblewhite et al., 2018). More distinct differences in density-dependent responses between migrant and resident large herbivores may be evident when populations are closer to K under low predation (e.g., Mysterud et al., 2011).

Nevertheless, our study provides some unique insight into how seasonal foraging and grouping patterns may contribute to the directional patterns we report in switching migratory tactics with changes in density (Eggeman et al., 2016). During the recent decline of the Ya Ha Tinda elk population, we documented an overall switching rate of 15%/yr, with resident elk being more likely to switch to being a migrant when the elk population was high, whereas migrant elk were more likely to switch to being a resident when the population was low. We hypothesize that the propensity of a resident elk to switch its migratory tactic at a high population reflects the effect of a lower quality summer range (Hebblewhite and Merrill, 2009) combined with density-dependent constraints on winter foraging, which is consistent with condition-dependent switching due to lower body reserves in spring when the decision to migrate occurs. An increase in the propensity for resident elk to migrate when the previous summer’s precipitation was low (Eggeman et al., 2016) supports this hypothesis. In contrast, the increased propensity of a migrant to remain as a resident at low density may stem not only from relaxed density-dependent competition (Mysterud et al., 2011, this study), but may be socially facilitated especially in open areas where herbivores tend to aggregate under high predation (Couzin et al., 2011; Lesmerises et al., 2018). Individuals in large groups would increase their conspecific interactions and social bonding without the same scramble competition for resources. In particular, if a migrant loses a calf to predation during the previous summer, the security of the human refuge and social ties on the winter range may foster the shift in migratory tactic (Switzer, 1997; Hoover, 2003; Gehr et al., 2020).

Our study fills a gap in the partial migration literature by pointing to key behaviors that may help explain the dynamics of different migratory tactics leading to the maintenance or loss of partial migration (Berg et al., 2019). Nevertheless, it falls short of linking these behaviors to reproductive outcomes. Previous reproduction within the context of spatial patterns of seasonal predation in conjunction with weather conditions may be key factors if the propensity to migrate is condition-dependent. Even where conservation efforts can prevent major habitat loss and barriers to migration, we expect migration patterns of ungulates to be dynamic in the face of large carnivore recovery and environmental change, and a focus on animal behavior may serve as a first indicator of what is to come (Sih, 2013; Greggor et al., 2016).

Data Availability Statement

The datasets generated for this study will not be made publicly available. Global Positioning System (GPS) and VHF telemetry data are available on the Movebank Data Repository: Hebblewhite et al. (2020).

Ethics Statement

The animal study was reviewed and approved by Animal Care Protocols: University of Alberta 353212, 611812, 000624, University of Montana AUP 004-16, Parks Canada BAN-2014-16756.

Author Contributions

JK, JP, MT, HM, and EM analyzed data and wrote sections of manuscript drafts. JB, HB, SE, and MH collected data and provided input to data analysis and manuscript writing. All authors approved this publication.

Funding

This work was supported by funding from Banff National Park, Alberta Fish and Wildlife, National Science Foundation LTREB grant 1556248, The Natural Sciences and Engineering Research Council of Canada (NSERC Discovery to EM), Montana Institute of Ecosystems (supported by National Science Foundation EPSCoRR RII Track 1 award), Universities of Alberta, and Montana, Alberta Conservation Association, Shikar Safari Club, International, Safari Club International, Safari Club International-Northern Alberta Chapter, Rocky Mountain Elk Foundation, Alberta Fish and Game and Hunting for Tomorrow for Ministers Special License funds.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank to M. J. Nyland, R. J. Smith, C. White, J. Whittington, B. Fyten, D. Gummer, B. Hunt, T. Shury, K. Heuer, T. Hurd, many other Parks Canada and Ya Ha Tinda ranch staff, and E. Bruns, J. Allen, and A. Hubbs from Alberta Fish and Wildlife, and for supporting this long-term project. We acknowledge the field help of C. Suthmeier, T. Weeks, M. Hessami, S. Hazenberg, S. Wales, C. Intering, C. Roberge, A. Bohm, S. L. Glines, P. Smolko, A. Deedy, and many volunteers.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00269/full#supplementary-material

References