DOI: 10.5061/dryad.2bvq83c0d
title: Data from: The primacy of density-mediated indirect effects in a community of wolves, elk, and aspen

description: 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.
description: <strong>Aspen data</strong>
   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.
   <strong>Elk
 data</strong>   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).   
 <strong>Wolf data</strong>
   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.  
 <strong>Weather data</strong>
   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.   
 <strong>Spatiotemporal variation in wolf predation
 risk</strong>  
 <em>Winter wolf spatial density</em>
   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.   <em>Kill
 Spatial Density</em>  
 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.  
 <em>Topography and vegetation
 openness</em>   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.
description: # Data from: The primacy of density-mediated indirect effects in a
 community of wolves, elk, and aspen
 [https://doi.org/10.5061/dryad.2bvq83c0d](https://doi.org/10.5061/dryad.2bvq83c0d) The file "BriceEtAl2024_Data.xlsx" includes all data needed to replicate models of aspen browse and recruitment probability as part of a test on density- and trait-mediated trophic cascades in Yellowstone National Park. The file "BriceEtAl2024_SupportingTables.xlsx" contains six supplemental tables from our analysis. ## File Structure ### BriceEtAl2024\_Data.xlsx This file has data on 26,012 individual young aspen collected between 1999-2019. Data were not collected in 2000 or 2015. The file has 3 tabs: one for data on aspen browse, and two for data on aspen recruitment. For the browsing data, each row represents an individual aspen at a specific site (i.e., "Plot", Column A), and the columns have data on whether that individual aspen was browsed by elk, and what it's height was. There are 10 columns for 10 separate measures of wolf predation risk of elk (described below and in methods), a column for annual elk abundance and annual elk density (abundance divided by 995 km2, i.e., study area), and columns for summer precipitation, winter snow water equivalent (SWE), and whether the site was within the 1988 Yellowstone fire boundary. Finally, there are three additional columns for height that were used in modeling: the first is a version of height that was scaled to have a mean of 0 and standard deviation of 1, and the second two are values of scaled height that are needed to model height as a spline. In the tabs for aspen recruitment, each row represents a single site in a single year, with the total number of aspen individuals that were counted at that site that year, and how many of them have "recruited" into the overstory, meaning they have a height taller than 120-cm. These tabs also have data for snow water equivalent, annual elk density, precipitation, and a single measure of wolf predation risk of elk (described below). Below, we describe each column in the dataframes. Note that "winter" refers to Nov 1st – Apr 30th.  #### Browse All data needed to reproduce models of browse probability in Brice et al. 2024. Columns are as follows: 1. Plot (Column A): ID number for our 113 aspen plots in northern Yellowstone National Park. 2. Year (Column B): the year that the data was collected. 3. Browse (Column C): binary variable for whether an individual aspen was browsed by elk the previous winter. A value of 0 indicates that the aspen was not browsed, and a value of 1 indicates that the aspen was browsed. 4. Height (Column D): height (in centimeters) of the tallest stem (i.e., leader stem) of individual aspen. 5. Wolf-Density-Aggregate (Column E): long-term average (i.e., 1999-2019) winter wolf density (wolves/30-m2) estimated for each site. For example, the value of 0.0000903 for plot 1 in 1999 means that, within a 30-m2 area, there was an estimated average of 0.0000903 wolves across all years of the study at plot 1. 6. Wolf-Density-Annual (Column F): winter wolf density (wolves/30-m2) estimated at each site for the corresponding year (Col B). For example, the value of 0.0000429 for plot 1 in 1999 means that, within a 30-m2 area, there was an estimated 0.0000429 wolves in 1999 at plot 1. 7. All-Elk-Kills-Aggregate (Column G): long-term average (i.e., 1999-2019) density (kills/30-m2) of all elk killed by wolves in winter at each site. For example, the value of 2.22E-09 for plot 1 in 1999 means that, within a 30-m2 area, there was an estimated average of 2.22E-09 elk killed by wolves across all years of the study at plot 1. 8. All-Elk-Kills-Annual (Column H): density (kills/30-m2) of all elk killed by wolves in winter at each site for the corresponding year (Col B). For example, the value of 2.67E-09 for plot 1 in 1999 means that, within a 30-m2 area, there was an estimated of 2.67E-09 elk killed by wolves in 1999 at plot 1. 9. Male-Elk-Kills-Aggregate (Column I): long-term average (i.e., 1999-2019) density (kills/30-m2) of adult and yearling male elk killed by wolves in winter at each site. This column can be interpreted the same way as Column G, but for just male elk instead of all elk. 10. Male-Elk-Kills-Annual (Column J): density (kills/30-m2) of adult and yearling male elk killed by wolves in winter at each site for the corresponding year (Col B).This column can be interpreted the same way as Column H, but for just male elk instead of all elk. 11. Female-Elk-Kills-Aggregate (Column K): long-term average (i.e., 1999-2019) density (kills/30-m2) of adult female elk and elk calves killed by wolves in winter at each site. This column can be interpreted the same way as Column G, but for just female elk and elk calves instead of all elk. 12. Female-Elk-Kills-Annual (Column L): density (kills/30-m2) of adult female elk and elk calves killed by wolves in winter at each site for the corresponding year (Col B). This column can be interpreted the same way as Column H, but for just female elk and elk calves instead of all elk. 13. Openness (Column M): proportion of aspen plot that was not tree cover, estimated for 30-m2 cells. A value of 0 indicates the landscape is completely covered by trees (i.e., all trees, 0% open), and a value of 1 indicates the landscape is completely open (i.e., no trees, 100% open) 14. Smoothness (Column N): how smooth/flat the topography is. This value is a relative scale, with no units. A value of 0 indicates very rough terrain, and a value of 1 indicates very smooth/flat terrain.  15. Elk-AnnualCount (Column O): annual counts of the elk population for each year. 16. Elk_Density (Column P): annual counts of the elk population divided by the study area to approximate elk density (elk/km2). For example, a value of 13.49 in 1999 at all stands indicates that there were an estimated 13.49 elk per km2 in 1999 17. Precipitation_cm (Column Q): sum of daily precipitation (cm; sum of all forms converted to water-equivalent) from Apr 1st – July 31st at each stand each year. For example, a value of 17.04 at plot 1 in 1999 indicates that there was an estimated 17.04 cm of rain at plot 1 during the summer of 1999. 18. SWE_tons (Column R): total winter snow water equivalent (tons/m2) at each stand each year. For example, a value of 14.49 at plot 1 in 1999 indicates that there was an estimated 14.49 tons of snow per m2 that accumulated during the winter before aspen was measured in 1999 at plot 1.  19. Fire (Column S): a binary variable indicating whether a plot was within the perimeter of the 1988 Yellowstone fires. A value of 1 indicates that the plot was within the fire perimeter, and a value of 0 indicates the site was outside of the fire perimeter. 20. Height-Scaled (Column T): aspen heights (Col D) scaled to have a mean of 0 and standard deviation of 1. This column was used to create the values needed for the height spline (Cols U & V) 21. Height-Spline1 (Column U): values for the first height spline, with a knot at 120-cm height. Aspen stems that are greater than or equal to 120-cm will have a value of 0.275. This column was used (along with Col V) to model aspen height as a non-linear predictor of aspen browse. 22. Height-Spline2 (Column V): values for the second height spline, with a knot at 120-cm height. Aspen stems greater than 120-cm will have a value of Height-Scaled - 0.275. Aspen stems that are less than 120-cm will have a value of 0. #### Recruitment There are two tabs of data for the recruitment models: (1) Annual wolf density and (2) Annual male elk kills. The recruitment data was created by counting the number of stems at each site each year that were taller than or equal to the previously estimated peak browse threshold, which was 120-cm in the annual wolf density model, and 123-cm in the annual male elk kill model. Note that the two tabs have different numbers of rows; this difference is because outliers in the risk variable were removed prior to estimating the number of recruited stems (i.e., some plot-year combinations missing from each tab). The columns for both tabs are all the same except for Column H, which is the risk variable. Descriptions of each column are as follows: 1. Plot (Column A): ID number for our 113 aspen plots in northern Yellowstone National Park. 2. Year (Column B): the year that the data was collected. 3. N_stems (Column C): total number of aspen stems counted at each site each year (regardless of height). For example, at plot 1, there were 37 stems measured in 1999 (row 2, col C), and 42 stems measured in 2001 (row 3, col C). 4. N_Recruits (Column D): total number of aspen stems taller than the peak browse threshold (i.e., 120-cm for annual wolf density model, 123-cm for annual male elk kill density model) at each site each year. For example, at plot 1, there were 0 stems taller than 120-cm in 1999 (row 2, col C, "Annual wolf density" tab), and 15 stems taller than 120-cm in 2010 (row 9, col C). 5. SWE_tons (Column E): total winter (Nov 1st – Apr 30th) snow water equivalent (tons/m2) at each stand each year. 6. Elk-Density (Column F): annual counts of the elk population divided by the study area to approximate elk density (elk/km2). 7. Precip_cm (Column G): sum of daily precipitation (cm; sum of all forms converted to water-equivalent) from Apr 1st – July 31st at each stand each year 8. Risk variable (Column H): 1. Wolf-Density-Annual: wolf density (wolves/30-m2) estimated at each site for the corresponding year 2. Male-Elk-Kills-Annual: density (kills/30-m2) of adult and yearling male elk killed by wolves in winter at each site for the corresponding year ### BriceEtAl2024\_SupportingTables.xlsx This excel file contains six tabs of supplementary results, each with it's own table. Each tab contains a text box with a description of the data. The tabs are as follows: 1. Table S1_Random Effects: when modeling aspen browse probability, we used generalized linear mixed models that contained random effects for plot and year. We also tested various combinations of random slopes to determine what random effect structure was the best fit for the model. This table shows the results from the model selection process.  1. Model (Column A): explains the random effect structure being tested. We tested various structures for each predation risk metric 2. N (Column B): Sample size for each model 3. Constant (Column C): the number of constants (i.e., intercepts) in the model. There was one constant in each model. 4. Rand. Intercept (Column D): the number of random intercepts included in the random effect structure. Two random intercepts (one for plot, one for year) were included in each model. 5. Rand. Slopes (Column E): the number of random slopes (and associated covariances) being tested in each model. A value of 0 indicates that there were no random slopes included in the model. In a model where we test the effect of including a random slope for elk density that varies by plot, we would have a random slope value of 2. This number includes the variance of the random slope for elk, as well as the covariance between elk and plot. In a model where we we test the effect of including a random slope for elk density that varies by plot AND a random slope for predation risk that varies by plot, we would have a random slope value of 5. This number includes the variance of the random slope for elk density, the variance of the random slope of risk, the covariance between elk density and plot, the covariance between risk and plot, and the covariance between elk density and risk. 6. Covariate (Column F): the number of covariates, or fixed effects, in the model. Each model included 5 fixed effects. 7. K (Column G): the total number of parameters in the model, which is the sum of columns C-F. 8. LL (Column H): the log likelihood of the model. 9. AICc (Column I): the AICc value for the model. The model with the lowest value is the best fit for the data. 10. ∆AICc (Column J): the difference in AICc between each model and the best model in the set. A value of 0 indicates that model is the best fit for the data. 11. Wi (Column K): the Akaike weight of the model. The Akaike weight of model *i* can be interpreted as the probability that model *i* is the best model given the collection of models considered. 2. TableS2_Height Spline: when modeling aspen browse probability, we included aspen height as a predictor. However, the effect of height on browse is nonlinear, with browse increasing with height up to a certain point, and then decreasing with height. To determine at which height browse switches from increasing to decreasing, we modeled height as a piecewise linear spline. This table shows the model selection results from this process. Each row is for a model with a different height threshold being tested. Cells highlighted in green are all within 2 AICc of the best model, which is indicated with bold font. The columns are similar to Table S1: 1. Knot (Column A): indicates what height was used as the point at which browsing switches from increasing with height to decreasing. 2. N (Column B): Sample size for each model. 3. Constant (Column C): the number of constants (i.e., intercepts) in the model. There was one constant in each model. 4. Rand. Intercept (Column D): the number of random intercepts included in the random effect structure. Two random intercepts (one for plot, one for year) were included in each model. 5. Rand. Slope (Column E): the number of random slopes (and associated covariances) being tested in each model. See item 1.5 above for further detail. 6. Covariate (Column F): the number of covariates, or fixed effects, in the model. Each model included 5 fixed effects. 7. K (Column G): the total number of parameters in the model, which is the sum of columns C-F. 8. LL (Column H): the log likelihood of the model 9. AICc (Column I): the AICc value for the model. The model with the lowest value is the best fit for the data 10. ∆AICc (Column J): the difference in AICc between each model and the best model in the set. A value of 0 indicates that model is the best fit for the data. 11. Wi (Column K): the Akaike weight of the model. 12. RE Structure (Column L): the random effect structure included in the model. 3. TableS3_Risk Shape: when modeling aspen browse probability, we wanted to know whether the effect of predation risk on browsing was linear, or if there was a threshold of risk, wherein there was little effect of risk up to a point, after which browsing decreased with risk. To test for a threshold, we compared models with a linear effect of predation risk to models with thresholds at low, moderate, or high risk, which were represented with piecewise linear splines. This table shows the model selection results from this process. Each row is for a model with a different risk threshold being tested. The columns are similar to Table S1: 1. Model (Column A): indicates the functional form of risk (i.e., linear, spline) being tested. 2. N (Column B): Sample size for each model 3. Covariate (Column C): the number of covariates, or fixed effects, in the model. Models with a linear form of risk have 5 fixed effects, while models testing for a threshold have two covariates for risk, and therefore have 6 fixed effects. 4. K (Column D): the total number of parameters in the model. 5. LL (Column E): the log likelihood of the model. 6. AICc (Column F): the AICc value for the model. The model with the lowest value is the best fit for the data. 7. ∆AICc (Column G): the difference in AICc between each model and the best model in the set. 8. Wi (Column H): the Akaike weight of the model. 4. TableS4_Elk Shape: when modeling aspen browse probability, we wanted to know whether the effect of elk density on browsing was linear, or if there was a threshold of elk density, below which there was a constant effect of elk density. To test for a threshold, we compared models with a linear effect of elk density to models with a log effect of density, and thresholds at 2-6 elk/km2, modeled with piecewise linear splines. This table shows the model selection results from this process. Each row is for a model with a different elk density threshold being tested. The columns are similar to Table S1: 1. Model (Column A): indicates the functional form of elk density (i.e., linear, log, spline) being tested. 2. N (Column B): Sample size for each model. 3. Rand. Intercept (Column C): the number of random intercepts included in the random effect structure. Two random intercepts (one for plot, one for year) were included in each model. 4. Rand. Slope (Column D): the number of random slopes (and associated covariances) being tested in each model. 5. Covariate (Column E): the number of covariates, or fixed effects, in the model. Models with a linear form of elk density have 5 fixed effects, while models testing for a threshold have two covariates for elk density, and therefore have 6 fixed effects. 6. K (Column F): the total number of parameters in the model, which is the sum of columns C-F. 7. LL (Column G): the log likelihood of the model. 8. AICc (Column H): the AICc value for the model. The model with the lowest value is the best fit for the data. 9. ∆AICc (Column I): the difference in AICc between each model and the best model in the set. 10. Wi (Column J): the Akaike weight of the model. 11. RE Structure (Column K): the random effect structure included in the model. 5. TableS5_Browse Models: coefficient estimates from the final models of browse probability as a function of height, SWE, elk density, and predation risk. Height was represented as a piecewise linear spline, with separate effects for stems shorter than ~120-cm than those taller than ~120-cm. The heights used for the spline are specific to each risk variable (see Table S4).  1. Risk metric (Column A): the predation risk metric included in the model. 2. Parameter (Column B): the parameters (i.e., predictors) included in the model. 3. β (Column C): the beta coefficient estimated for each parameter in the model. 4. 95% Confidence Interval (Columns D & E): the lower (col D) and upper (col E) limits of the 95% confidence interval for the parameter estimate. 5. SE (Column F): the standard error estimated for the parameter. 6. P (Column G): the p-value associated with the parameter. 7. RE Structure (Column H): the random effects structure included in the model. 6. TableS6_Density-Risk Interaction: coefficient estimates from models of browse probability as a function of height, SWE, elk density, predation risk, and an interaction between elk density and risk. Height was represented as a piecewise linear spline, with separate effects for stems shorter than ~120-cm than those taller than ~120-cm. The heights used for the spline are specific to each risk variable (see Table S4).  1. Risk metric (Column A): the predation risk metric included in the model. 2. Parameter (Column B): the parameters (i.e., predictors) included in the model. 3. β (Column C): the beta coefficient estimated for each parameter in the model. 4. 95% Confidence Interval (Columns D & E): the lower (col D) and upper (col E) limits of the 95% confidence interval for the parameter estimate. 5. SE (Column F): the standard error estimated for the parameter. 6. P (Column G): the p-value associated with the parameter. 7. RE Structure (Column H): the random effects structure included in the model.