The primacy of density-mediated indirect effects in a community of wolves, elk, and aspen

The removal or addition of a predator in an ecosystem can trigger a trophic cascade, whereby the predator indirectly influences plants and/or abiotic processes via direct effects on its herbivore prey. A trophic cascade can operate through a density-mediated indirect effect (DMIE), where the predator reduces herbivore density via predation, and/or through a trait-mediated indirect effect (TMIE), where the predator induces an herbivore trait response that modifies the herbivore’s effect on plants. Manipulative experiments suggest that TMIEs are an equivalent or more important driver of trophic cascades than are DMIEs. Whether this applies generally in nature is uncertain because few studies have directly compared the magnitudes of trait- and density-mediated indirect effects on natural unmanipulated field patterns. A TMIE is often invoked to explain the textbook trophic cascade involving wolves (Canis lupus), elk (Cervus canadensis), and aspen (Populus tremuloides) in northern Yellowstone National Park. This hypothesis posits that wolves indirectly increase recruitment of young aspen into the overstory primarily through reduced elk browsing in response to spatial variation in wolf predation risk rather than through reduced elk population density. To test this hypothesis, we compared the effects of spatiotemporal variation in wolf predation risk and temporal variation in elk population density on unmanipulated patterns of browsing and recruitment of young aspen across 113 aspen stands over a 21-year period (1999-2019) in northern Yellowstone National Park. Only two of ten indices of wolf predation risk had statistically meaningful effects on browsing and recruitment of young aspen, and these effects were 8-20 times weaker than the effect of elk density. To the extent that temporal variation in elk density was attributable to wolf predation, our results suggest that the wolf-elk-aspen trophic cascade was primarily density-mediated rather than trait-mediated. This aligns with the alternative hypothesis that wolves and other actively hunting predators with broad habitat domains cause DMIEs to dominate whenever prey, such as elk, also have a broad habitat domain. For at least this type of predator-prey community, our study suggests that risk-induced trait responses can be abstracted or ignored while still achieving an accurate understanding of trophic cascades.

Methods
Aspen data

Beginning in 1999, we measured young aspen height and browsing at 113 stands selected in a stratified random sample reflecting high and low wolf use areas (see Brice, Larsen, and MacNulty 2022 for details). All stands were selected from aerial photographs taken after the 1988 fires; as such, selected aspen stands were those whose overstory at least partially survived the 1988 fires. Each stand contained a 20-x-1 m belt transect, and we surveyed all young aspen (≤ 600 cm tall, “stems”) within the transect (“plot”) at the end of the growing season (late July – September) each year. For each stem, we measured height of the leader (i.e., tallest stem), and whether the leader was browsed the previous winter. The number of stands sampled each year varied from 61 – 113 (μ = 97.3, σ = 18.3), and the number of plots with stems each year varied from 55 – 108 (μ = 84.9, σ = 17.2). Sampling occurred annually from 1999-2019, excluding 2000 and 2015, resulting in 26,012 stem-level observations over 19 years.

Elk data

Aerial winter surveys of elk were conducted annually using 3-4 fixed wing aircraft simultaneously flying non-overlapping areas between Dec – Mar (see Lemke, Mack, and Houston 1998). In years with no survey (i.e., 2006, 2014), elk counts were interpolated with a state-space model and corrected for the effects of elk group size on sightability (Tallian et al. 2017; B. J. Smith and MacNulty 2023). We divided annual counts of elk within the Park by the study area (995 km2) to calculate annual elk density (number of individuals per km2). 

Wolf data

Since 1995, the Yellowstone Wolf Project has studied wolves for two 30-day periods each winter: (1) mid-Nov to mid-Dec (early winter) and (2) the month of March (late winter). Each winter, 20-30 wolves (~35-40% of population) were captured and fitted with VHF and GPS collars (D. W. Smith et al. 2004). All wolves were captured and handled following protocols in accordance with guidelines from the American Society of Mammalogists (Sikes 2016) and approved by the National Park Service Institutional Animal Care and Use Committee (IACUC permit IMR_YELL_Smith_wolves_2012). All wolf packs in northern YNP had at least one collared wolf each year. Locations from both VHF and GPS collars were recorded approximately daily during early and late winter periods, and weekly outside of these periods. GPS collars also recorded hourly locations during each 30-day winter study, and at variable times otherwise. During winter study, ground and aerial crews searched for wolf kills by tracking collared wolves and investigating clusters of locations.

Weather data

We obtained data on SWE at each aspen stand from Daymet, which produced daily gridded estimates of weather parameters from meteorological observations at a 1-km2 resolution (Thornton et al. 2020). We calculated total winter SWE (tons/m2) by summing daily estimates from Nov 1st – Apr 30th at each stand each year. We also obtained data on spring precipitation, which we estimated as the sum of daily precipitation (cm; sum of all forms converted to water-equivalent) from Apr 1st – July 31st, again obtained from Daymet for each stand each year. 

Spatiotemporal variation in wolf predation risk

Winter wolf spatial density

We used VHF and GPS locations of wolves in the study area to calculate winter (Nov 1 – Apr 30) wolf density each year and across years. We restricted the data to wolves with at least 30 days of observations, which proved to be highly correlated with the full 6-months of locations (Pearson’s r = 0.99). Additionally, we only used wolves with at least 10 locations per winter, the minimum number of locations needed for the models to converge. After restricting the data, there were 142,087 total locations and 777 unique wolf-year combinations (“wolf-years”) from 1999 – 2019, with wolf-years spanning 30 – 181 days (median = 152 days) and containing 10 – 4,194 locations (median = 42).

To estimate the spatial densities of wolves, we used the locations to fit individual continuous time movement models (CTMM) to each wolf-year using the ctmm package (Calabrese, Fleming, and Gurarie 2016) in R (V1.2.5019, R Core Team 2018). We used the Ornstein-Uhlenbeck Foraging (OUF) anisotropic process for each wolf, which accounts for correlated velocities and restricted space use (C. H. Fleming et al. 2014). Once each wolf had its own CTMM, we calculated an autocorrelated kernel density estimate (AKDE) at a 30-m2 resolution for each wolf-year. If there were multiple collared wolves within a pack, we averaged their AKDEs and divided by the sum of all values to ensure that the AKDE summed to one and could be interpreted as a probability density.

Once we had a single AKDE for each pack each winter, we weighted each pack-specific AKDE by the corresponding number of wolves in each pack (lone wolves unweighted), and then summed the densities of all packs and lone wolves each winter, resulting in a single wolf AKDE each year. Finally, we created a long-term average measure of wolf density by taking the mean of all annual AKDEs. We intersected all spatial layers of risk with the aspen stand locations to derive stand-specific estimates of risk.

Kill Spatial Density

We used positional data of wolf-killed elk to calculate a kernel density estimate (KDE) of elk kills each winter using the sp.kde function from the spatialEco package in R (Evans, Murphy, and Ram 2021). We used a bandwidth of 3 km per the methods of Kohl et al. (2018) and Fortin et al. (2005), and a resolution of 30-m2 (Kauffman et al. 2007; Kohl et al. 2018). We distinguished kills by sex, creating annual KDEs with all kills (N = 2448, Annual range = 61-193), adult male elk and male yearlings (N = 729, Range = 17-69), and adult female elk and calves (N = 1430, Range = 28-125) for each winter. As with wolf density, we also calculated long-term averages of kill density using kills across all years for the three categories.

Topography and vegetation openness

We extracted land-cover type using the Rangeland Analysis Platform (Allred et al. 2021), and calculated openness for each year as the proportion of each 30-m2 cell that was not tree cover. To calculate smoothness, we used a 30-m2 digital elevation model (DEM) and the terrain function from the raster package in R (Hijmans and van Etten 2014), which produced a map of roughness. Roughness was defined as the difference between the maximum and minimum elevation of a cell and its surrounding 8 cells. We converted roughness to smoothness by scaling it from 0-1 and subtracting from 1.