1. Ecology
Download icon

Linking spatial patterns of terrestrial herbivore community structure to trophic interactions

  1. Jakub Witold Bubnicki  Is a corresponding author
  2. Marcin Churski
  3. Krzysztof Schmidt
  4. Tom A Diserens
  5. Dries PJ Kuijper
  1. Mammal Research Institute, Polish Academy of Sciences, Poland
Research Article
  • Cited 0
  • Views 833
  • Annotations
Cite this article as: eLife 2019;8:e44937 doi: 10.7554/eLife.44937

Abstract

Large herbivores influence ecosystem functioning via their effects on vegetation at different spatial scales. It is often overlooked that the spatial distribution of large herbivores results from their responses to interacting top-down and bottom-up ecological gradients that create landscape-scale variation in the structure of the entire community. We studied the complexity of these cascading interactions using high-resolution camera trapping and remote sensing data in the best-preserved European lowland forest, Białowieża Forest, Poland. We showed that the variation in spatial distribution of an entire community of large herbivores is explained by species-specific responses to both environmental bottom-up and biotic top-down factors in combination with human-induced (cascading) effects. We decomposed the spatial variation in herbivore community structure and identified functionally distinct landscape-scale herbivory regimes (‘herbiscapes’), which are predicted to occur in a variety of ecosystems and could be an important mechanism creating spatial variation in herbivory maintaining vegetation heterogeneity.

https://doi.org/10.7554/eLife.44937.001

eLife digest

In almost every ecosystem on Earth, communities of herbivores are kept in check by both predators and the availability of the plants they eat. As herbivores move in response to these pressures, they shape local plant communities and impact vegetation across entire landscapes. Yet the role of large plant-eating mammals in structuring ecosystems is often overlooked. Indeed, most research on this topic has looked at African ecosystems, like open savannahs, and fewer researchers have studied temperate forests like those found across Europe, Asia and North America.

Bubnicki et al. have now examined factors influencing the distribution of five large herbivore species and resulting plant communities in Białowieża Forest in eastern Poland, the best-preserved European lowland forest. Their method involved measuring the cascading interactions of plants and animals in the forest using cameras set at nearly 900 locations, satellite images and other remote sensing technologies, and on-the-ground surveys. Added to this were patterns of human activity inferred from the available data for the study area. This approach allowed Bubnicki et al. to explore how humans are influencing the forest ecosystem, too.

The analysis revealed that humans are the main factor influencing the movements of carnivorous predators in Białowieża Forest, but not the herbivores directly. Wolves and lynxes avoided areas heavily used by humans whereas large herbivores responded primarily to different environmental factors. Wild boar and bison are influenced by the availability of plant food and preferred habitat for foraging; moose and roe deer by the features of the landscape, like elevation or openness. The red deer was the only large herbivore species whose distribution was strongly linked to that of its main predator, the wolf.

From this, Bubnicki et al. identified distinct areas in the forest which have emerged from the interactions at play, describing these areas as ‘herbiscapes’ for the herbivores that shaped them. These findings provide new understanding of the complex ecological processes shaping the Białowieża Forest and serve as a model to help understand other ecosystems around the world. The knowledge will also contribute to the ongoing management and conservation of this UNESCO World Heritage Area.

https://doi.org/10.7554/eLife.44937.002

Introduction

Spatial patterns in species distribution, abundance and community composition are manifestations of underlying ecological mechanisms operating at a range of spatial scales (Levin, 1992). These patterns emerge from dynamic interactions between environmental bottom-up, biotic top-down and biotic parallel factors in combination with stochastic effects. Although the role of bottom-up factors in shaping species distributions has been intensively studied in recent decades (Elith and Leathwick, 2009; Guisan and Thuiller, 2005), we still know little about the importance of the biotic factors underlying most spatial patterns. Especially, the role of species interactions, within and across trophic levels, including those involving humans, remain largely unexplored (Darimont et al., 2015; Schmitz et al., 2017; Wiens, 2011; Worm and Paine, 2016). It has recently been proposed that community ecology should be ‘rediscovered’ as an integrative study of species interactions and spatial distributions (Schmitz et al., 2017), while accounting for direct and indirect anthropogenic effects on species distributions and behavior (Berger, 2007; Worm and Paine, 2016).

Large mammalian herbivores influence terrestrial ecosystem structure and functioning (Gordon et al., 2004; Hobbs, 1996; Schmitz, 2008) via their direct effects on vegetation structure (Charles-Dominique et al., 2016; Churski et al., 2017; Didion et al., 2009; Hempson et al., 2015; Kuijper et al., 2010a) and indirect effects on nutrient cycling (Murray et al., 2013). In this way, herbivory influences vegetation at large spatial scales, from the local landscape up to the biome level (Moncrieff et al., 2016; Woodward et al., 2004), and can lead to herbivory-mediated cascading effects on other trophic levels (Gordon et al., 2004; Palmer et al., 2015; Schmitz, 2008). There is often strong spatial variation in herbivory impact, resulting from the space use of different functional groups of herbivores. This spatial variation is driven primarily by the interactive effects of abiotic factors, disturbances, forage quality and quantity in combination with life-history traits, such as herbivore body mass (Anderson et al., 2016; Cromsigt et al., 2009; Hempson et al., 2015; Hopcraft et al., 2010; Ogutu et al., 2010). However, the actual distribution of many herbivores often differs from the expected distribution derived purely from interactions with bottom-up factors. This discrepancy results from herbivores also responding to landscape gradients induced by biotic top-down interactions (Anderson et al., 2010; Hopcraft et al., 2010; Kauffman et al., 2007). In effect, this landscape of interacting ecological gradients, both bottom-up and top-down, creates spatial heterogeneity in the availability and suitability of habitats for different large herbivore species within a community (Cromsigt et al., 2009; Fryxell, 1991; Hopcraft et al., 2010). Thus, to assess the ecosystem-level impact of large herbivore communities requires full understanding of the factors driving spatial heterogeneity in their community structure across a landscape (Gordon et al., 2004; Weisberg and Bugmann, 2003).

Recently, there has been much attention given to the role of large carnivores in structuring ecosystems via their effects on herbivore communities (Estes et al., 2011; Ripple et al., 2014; Terborgh et al., 2006). In addition to their density-mediated effects (i.e. impact on prey population size), behaviorally mediated effects (i.e. impact on prey behavior) on prey species are a crucial mechanism explaining the trophic cascades driven by large carnivores (sensu Ripple et al., 2016). Prey species react to the presence of large carnivores by adjusting their spatio-temporal patterns of landscape use (Creel et al., 2005; Kohl et al., 2018; Laundré et al., 2001; Valeix et al., 2009). These spatial interactions between trophic levels are usually context-dependent and are shaped by the biophysical characteristics of a landscape (Kauffman et al., 2007; Schmitz et al., 2017; Valeix et al., 2009). Many studies have addressed how carnivores affect their prey species, but these have generally used a single carnivore - single prey species approach, whereas many ecosystems host multiple carnivore and multiple prey species. In such systems, different carnivore species can create contrasting risk effects (Creel et al., 2017; Preisser et al., 2007; Thaker et al., 2011). Moreover, in multi-species communities, some prey species perceive more risk than others from a carnivore species (Anderson et al., 2016; Laundré et al., 2001; Valeix et al., 2009).

This suggests that spatial distributions of predation-sensitive prey species may be mainly driven by carnivore top-down effects, whereas distributions of predation-insensitive prey or non-target species by gradients in resources availability (Hopcraft et al., 2010). When a community of prey species consists of ecologically or functionally similar species, changes in the abundance and distribution of one species may be buffered by another species (Ford et al., 2015; Rosenfeld, 2002). These so-called redundancy effects, can prevent apex predators from creating trophic cascading effects when taking the response of the entire herbivore community into account, despite them significantly impacting one or more prey species (Ford et al., 2015; Liu et al., 2016).

There is a growing awareness that including humans in community studies is critical for improving our understanding of ecosystem functioning (Darimont et al., 2015; Worm and Paine, 2016) and for predicting species distributions in increasingly anthropogenic environments. Due to the recovery of some large carnivore populations and expansion of human populations, carnivores are increasingly sharing landscapes with humans world-wide (Carter and Linnell, 2016; Chapron et al., 2014). Humans are also increasingly being considered a coherent part of complex trophic interaction chains (Darimont et al., 2015; Strong and Frank, 2010; Kuijper et al., 2016). The resulting complex, cascading interactions urgently need to be considered when studying the spatial distributions of herbivores, their effects within the landscape and the functional role large carnivores can play in landscapes that are becoming increasingly anthropogenic (Kuijper et al., 2016).

In this study, we investigated how the interactive effects of bottom-up and natural top-down factors (two large carnivore species and humans), determine the landscape distribution and community composition of five native ungulate species in Białowieża Forest (BF, Poland; Figure 1). BF is regarded to be one of the best preserved temperate European lowland forest systems and is inhabited by a natural community of large mammals (Jędrzejewska and Jędrzejewski, 1998). In addition, BF is also embedded within an anthropogenic landscape typical for many terrestrial systems. We hypothesized that spatial variation in the composition of the large herbivore community is explained by the interactive effects of species-specific responses to major environmental and risk gradients operating at the landscape level. We aimed to answer the following questions: 1) Do large carnivores have species-specific effects on the distributions of ungulates in our multiple predator-prey system?, 2) How does human activity mediate predator-prey interactions at the landscape scale?, 3) Does this lead to ecologically distinct herbivory regimes (sensu Hempson et al., 2015) with differential vegetation impact at the landscape scale? Using detailed data on species distributions (894 camera trap locations), landscape structure (high-resolution GIS and remote sensing data) and detailed woody vegetation surveys (385 study plots) along with a novel spatially-explicit hierarchical modelling approach we decomposed the spatial variation in herbivore community structure into ecologically distinct landscape-scale herbivory regimes. With data from a complex, multi-species and human-influenced system, we aimed to exemplify the functional role that large carnivores and their herbivorous prey can play in increasingly human-affected ecosystems.

Maps of the study area (left) and camera trapping sampling design (right).

On the study area map we marked the major settlements connected by public roads open for cars.

https://doi.org/10.7554/eLife.44937.003

Materials and methods

Study area

Request a detailed protocol

The study was carried out in Białowieża Forest (BF) in Poland (c. 580 km2; Figure 1). This harbors a natural assemblage of central European ungulate species, with red deer (Cervus elaphus) being the most abundant (6.0 individuals/km2), followed by wild boar (Sus scrofa; 5.4/km2) and roe deer (Capreolus capreolus; 2.0/km2); and European bison (Bison bonasus; 0.5/km2) and moose as the rarest ungulates (Alces alces; 0.08/km2) (Borowik et al., 2016). Two large carnivores occur in BF: the Eurasian lynx (Lynx lynx; c. 15 individuals) and wolf (Canis lupus; 4 packs of 7–12 individuals) (Jędrzejewski et al., 2002; Schmidt and Kuijper, 2015). Part of BF, the core of Białowieża National Park (BNP; c. 47 km2), has been strictly protected since 1921. Since this time, human activities such as hunting and forestry have been banned. In 1996 BNP was enlarged to c. 100 km2. Outside the national park, the forest is managed by the State Forest Holding, hence timber production takes place and ungulate hunting is allowed, but here also exists a network of nature reserves (c. 130 km2, see Figure 1). For a more detailed description of this area see Faliński and Falińska (1986) and Jędrzejewska and Jędrzejewski (1998). From a landscape perspective, it is important that the boundaries of BF ecosystem are well defined. From the west, BF is surrounded by agricultural fields, and from the north and south by a mosaic of agricultural and fragmented forest landscapes. In the east, a tall wire-fence along the border with Belarus (built in 1981) prevents movements by ungulates. These conditions create an opportunity to study the spatial distribution of a whole community of ungulates in a spatially restricted, complex ecosystem with natural predators present and varying levels of human management and impact.

Camera trapping

Request a detailed protocol

In contrast to telemetry, which allows the study of both predator and prey movement patterns by placing sensors on individual animals, we used point-based camera trapping to record and reconstruct the distributional patterns of large carnivores and ungulates. We sampled the continuous BF landscape with camera traps placed randomly with respect to species movements and hypothetical mechanisms driving their distributions (habitat structure, humans, predator-prey interactions). We used an intensive, large-scale and high-resolution camera trap network covering the entire study area to collect detailed, spatially-explicit information on species distributions. We argue that camera trapping is the most objective and efficient method for collecting this type of spatially-explicit community data. This kind of study would be inherently impossible to do with telemetry: to record the spatial patterns we were looking for, the entire populations of each species had to be observed simultaneously. Moreover, because of logistical and financial limitations and ethical issues (related to live-trapping of protected species in Europe) it would be practically impossible to obtain sufficient data for all the studied species of large herbivores and carnivores using telemetry. We used digital trail cameras (Ecotone SGN-5210A) triggered by passive infrared sensors with a detection angle of c. 35° and range of c. 20 m. After detection, with a time lag of c. 1 s, a photograph was taken and the camera recorded a 60 s video. When an animal stayed, this procedure was repeated without trigger delay. During low-light conditions, cameras switched to a stealth infrared mode. Cameras were attached to a tree at a height of c. 1 m at locations with a clear view of at least 20 m. Whenever possible, we randomly chose an acceptable place to mount a camera as close as possible to the locations given by coordinates pre-computed prior to the field work.

Ungulate survey

Request a detailed protocol

We quantified the spatial distribution of the ungulate community by using a spatially extensive, high-resolution network of camera traps. The data were collected over 2 years (May 2012 - May 2014) of intensive camera trapping in all seasons, except for days with the strongest winter conditions (snowing heavily or temperatures below −20°C). We conducted 34 trapping sessions, each lasting 14 days on average. During each session between 30 to 40 camera traps were pseudo-randomly deployed in different parts of the forest and within a minimum distance of 100 m from the nearest roads (both paved and unpaved), large clearings and settlements. Additionally, we kept a minimum distance of 100 m between all camera locations. The coordinates for all locations were pre-computed prior to the field work in QGIS software (QGIS Development Team, 2017). In total, we collected data at >1 k sites. However, because of logistical errors, camera failures, stealing of equipment and heavy snowing in winter (leading to blocked view and uninterpretable gaps in data) we had to exclude >100 sites (~10%) from further analysis. Finally, we used data collected at 894 sites, covering the whole BF landscape (Figure 1).

Carnivore survey

Request a detailed protocol

Large carnivores generally tend to have a higher detection probability at forest roads and trails (Cusack et al., 2015) because they often prefer to move along linear landscape structures (Zimmermann et al., 2014) and/or to increase the probability of encountering prey (Whittington et al., 2011). This specific space use resulted in a very low trapping rates of both carnivores (wolf and lynx) during the ungulate survey (Appendix 1—table 1), when camera traps were placed randomly in the forest. Hence, to better quantify the space use of the two large carnivores we ran an additional camera trapping survey in September-October 2015. We deployed 73 camera traps on the sides of forest roads across the whole landscape for one month (Figure 1). The core areas of wolf pack territories are related to the locations of breeding dens (Jędrzejewski et al., 2001). During the reproductive season (spring-summer) the spatial distribution of a wolf pack is restricted to their core area, whereas outside this period they regularly return to it (Jędrzejewski et al., 2001). In August-September, pups begin to travel with other pack members and move more widely through their territory, returning to the core of their territory on a regular basis (Jędrzejewski et al., 2001). Lynx have a similar, typical pattern of movements whereby females restrict their movements in May-July, while tending their kittens (Schmidt, 1998). Therefore, late summer-autumn is the best period to quantify the space use of both carnivores at the landscape scale.

The ungulate and carnivore surveys were conducted during different time periods, with different sampling intensities (May 2012 to May 2014 versus September to October 2015, respectively) and using different camera placement strategies. To check if our results were not spatio-temporally confounded because of these differences, we ran the same model for all ungulate species using only a subset of the camera trap data covering a 3-month period (August - October) matching the period of the carnivore survey as closely as possible (the models did not converge for data from only one or two months). The obtained results were similar to those based on the full dataset (Appendix 1: Figs. S36-S37); we thus chose to use the full dataset, which provided a larger sample size and better spatial coverage of the study area, both of which are needed for making robust inferences using complex hierarchical spatial models such as ours (see Statistical model section). Moreover, between-season variation in the trapping rate of ungulates was directly accounted for in the model. Lastly, all the studied ungulate species are non-migratory; in other words, there is no large (landscape)-scale seasonal movement of ungulates in BF (Jędrzejewska and Jędrzejewski, 1998; Kamler et al., 2008; Podgórski et al., 2013). The winter distribution of European bison is to a large extent driven by the location of supplementary feeding sites. This results in concentration of bison at these locations during the winter season.

Data processing

Request a detailed protocol

After downloading camera trap data, both ungulates and wolf datasets were organized and classified using TRAPPER software (Bubnicki et al., 2016). Species, sex, age and group size were determined for every recorded image or video containing an observation of focal species. We defined the independence interval between successive captures (i.e. event; see Meek et al., 2014) as five minutes.

Statistical model

General description

Request a detailed protocol

We developed a hierarchical multi-scale spatial model to quantify landscape use by large carnivores and ungulates. Our model was built upon the previous work of Royle (2004) and Royle et al. (2007), who described a class of Binomial-Poisson N-mixture models for spatially replicated counts collected during multiple (discrete) surveys at each site (Royle, 2004) and further applied these models in a spatially explicit context (Royle et al., 2007). Following the later work of Guillera-Arroita et al. (2011); Guillera-Arroita et al., 2012), we expanded this approach to a continuous case where counts are described using the Poisson instead of binomial distribution. This approach is more suitable for camera traps and unmarked populations (as in our case), as multiple detections of the same individuals are allowed, accounting within a single modeling framework for both false-negatives (‘imperfect detection’) and false-positives (‘double counts’).

Moreover, the Poisson distribution intensity parameter (γ in the next section), when multiplied by the number of (arbitrarily defined) sampling occasions, can be interpreted as the trapping rate, and also corrects for unequal effort amongst sampling locations (offset). Both the Binomial-Poisson N-mixture model and its continuous variations allow for estimating species abundances at monitored sites (Nj in the next section), assuming sampling locations are independent and that the observed system is closed (i.e. no changes in abundance at the site during repeated surveys). In other words, site abundances are treated as independent random variables distributed according to some mixing distribution, for example Poisson or Negative Binomial (Royle, 2004).

Usually, to ensure independence between sampling locations (i.e. no shared individuals) a distance larger than the average home range of focal species is preferred. This does not apply to our high-resolution camera trapping study as both the average distance between neighbouring sampling locations and size of landscape grid cells chosen for predictions (500 x 500 m) are much smaller than home and daily ranges of all ungulates and both carnivores. This specific sampling strategy was designed to capture and model the fine-scale, continuous variation in the use of the landscape by all studied species. Thus, following Royle (2004) we view the site-specific abundance as a random effect and relax the assumption of sampling site independence by interpreting Nj as the relative density that is the number of individuals using a given landscape grid cell during a sampling period rather than an absolute value of abundance. A similar interpretation has been used in many occupancy studies (Cusack et al., 2017; Efford and Dawson, 2012; Latif et al., 2016), where the probability of site occupancy has been interpreted as the probability of site use. We used the Negative Binomial distribution as the prior distribution for N to accommodate extra-Poisson variation not explained by the included covariates and spatial random effects (see Formal description below). Lastly, we assumed the system to be closed (i.e. no significant changes in ungulate and carnivore distributions and demographies) during the sampling period (2 years for ungulates and 1 month for carnivores). This is a reasonable assumption since the landscape surrounding BF creates natural boundaries (see description of the study area) and the hunting intensity is relatively low and constant over the years (Zbyryt et al., 2018).

Formal description

Request a detailed protocol

Let us consider a landscape divided into j=1,2,...,G grid cells of equal sizes and i=1,2,...,S sampling locations monitored with camera traps for di days each, where each grid cell can contain zero, one or multiple camera trap locations. This results in a multi-scale design where smaller subunits (camera trap sites) are nested within larger units (grid cells) (Kery and Royle, 2016). Counts, yi, are Poisson random variables given by

(1) yi|Nj[i]Poisson(Nj[i]γidi)

where the intensity parameter is a product of Nji, the number of individuals using a landscape grid cell j during a study period, γi, the expected detection (trapping) rate per sampling occasion (here 1 day) and di, the number of sampling occasions (days) during a survey at a camera trap location i. In the state-space formulation of our hierarchical model Nj is the Negative Binomial distributed latent variable

(2) NjNegBinλj,ϕ

with the parameter λj being the expected number of individuals using a landscape grid cell j and ϕ the dispersion parameter. Both parameters, λj and γi can depend on covariates describing for example environmental gradients (at different scales), biotic interactions and seasonal differences in species activity level. This variation can be modelled using standard generalized linear regression techniques using log-link functions

(3) logγi=xγi'βγ+ϵi
(4) logλj=xλj'βλ

where xγ' and xλ' are transposed rows of design matrices Xγ and Xλ, respectively, β are vectors of linear predictor coefficients and ϵi are identically and independently distributed iid camera trap measurement errors. It is necessary to note that, while the linear predictor of λj explains the variation in data arising from the ecological processes, the linear predictor of γi (detection/trapping rate) deals with both the ecological (e.g. habitat selection, movement, seasonal activity levels) and observational processes (e.g. detection issues, camera failures). In the latter case, informative and ecologically meaningful covariates are needed to disentangle these otherwise confounded sources of variation. In order to directly account for the potential spatial dependence between landscape grid cells (i.e. to capture all the spatial variation not explained by the included covariates) we introduced spatial random effects (“spatial residuals”) into the linear predictor of λj:

(5) log(λj)=xλjβλ+ωj

where ω=ω1,...,ωG is the realization of a Gaussian spatial process on a discrete (gridded) spatial domain. Specifically, we implemented a Restricted Spatial Regression model (RSR) (Hughes and Haran, 2013; Johnson et al., 2013), which is the restricted version of the intrinsic Conditional Autoregressive (iCAR) model. The RSR model was developed to solve the issue of confounding between the spatial process and the fixed-effects covariates in spatial regression models. The byproduct of this solution is its computational efficiency, as the RSR model is a reduced dimension model.

Specification

Request a detailed protocol

The model as described above was fitted to our camera trapping data of the wolf and lynx and each of the five ungulate species. To model the intensity of landscape use (λj in Equation 5), we overlaid a grid of 2303 cells 500 m per side (25 ha each) over the study area. Based on this grid we compiled a set of spatial (raster) covariates describing the ecologically relevant (sensu Elith and Leathwick, 2009) environmental and human-induced gradients (Figure 2). For GIS data processing, we used QGIS (QGIS Development QGIS Development Team, 2017) and GRASS GIS (Neteler et al., 2012) open source software. We standardized (scaled) all covariates, by subtracting the mean and dividing the result by the SD of the original variables. Finally, we ensured that the Pearson correlation coefficient for all pairs of included covariates was lower than 0.7 (Appendix 1—figure 1). For more details about all spatial (raster) covariates and the processing of GIS and remote sensing data see Appendix 1.

Maps of candidate landscape-scale covariates for the ungulate and carnivore models representing environmental and human-related landscape gradients.

All covariates were scaled and zero-centered prior to modelling. POIs – density of tourist infrastructure. When selecting covariates for the final models we ensured that the Pearson correlation value for all pairs of included covariates was lower than 0.7 (Appendix 1—figure 1).

https://doi.org/10.7554/eLife.44937.004

Specification - wolf and lynx model

Request a detailed protocol

Based on existing knowledge (Kuijper et al., 2015; Schmidt et al., 2009; Theuerkauf et al., 2003) both carnivore species in the Białowieża forest utilize the entire landscape, although with clear spatial patterns in the intensity of use. Previous work in this study area has shown that this is primarily determined by human related factors (Theuerkauf et al., 2003). In BF tourist traffic concentrates mainly within the central parts of the forest where roads (open for the public and cars) connect the three major settlements in the area, that is Hajnówka, Białowieża and Narewka (see Figure 1 and Appendix 1—figure 3). For the reasons above, the following raster layers were chosen as landscape covariates likely to influence carnivore space-use: distance to major settlements, distance to touristic trails, density of touristic infrastructure (POIs), density of protected areas (BNP and nature reserves) and elevation (Figure 2). Elevation was included as in this flat landscape the lowest, often swampy areas are the least accessible for humans and could therefore be preferred by the wolf and lynx. We interpreted the parameter λ as large carnivore space use intensity and included rasters with predicted values of λj for both species as covariates in all models for ungulates. We assumed that from a prey perspective, λj is proportional to the predator encounter rate, hence it quantifies potential risk as perceived by ungulates.

Specification - ungulate model

Request a detailed protocol

To explain the variation in the landscape-scale distribution of ungulates, we considered three environmental gradients primarily related to major biophysical properties of the forest environment that are known to affect space use of ungulates at multiple scales: percentage of landscape openness, tree stand canopy height and percentage share of coniferous species (Churski et al., 2017; Jedrzejewska et al., 1994; Kuijper et al., 2009; Kuijper et al., 2010a; Figure 2). For landscape openness and percentage share of coniferous species, we additionally included their quadratic effects, allowing for the existence of an optimum value for each variable (e.g. species preference for a mixed forest or for intermediate levels of canopy closure). The other covariates included were potential predation risk variables (space use of large carnivores) and two human-related landscape gradients, namely distance to all settlements and density of protected areas (BNP and nature reserves; Figure 2). Additionally, we computed the same set of environmental covariates but with a higher resolution (100 m) and used them to model variation in detection rate (γi in Equation 4) at camera trap sites. Here, we assumed that a part of this variation comes from a resource selection process operating at a scale smaller than the landscape unit we defined (see e.g. Johnson, 1980), influencing at-site detection rates and observed counts in the end. Another source of variation is species movement behavior and activity level, which can both change between seasons. To control for these temporal effects, we considered a quadratic function of temperature and snow cover as covariates for detection rate.

Implementation

Request a detailed protocol

The models were implemented within a Bayesian framework in Python using PyMC 2.3.6 software (Patil et al., 2010). We used Markov chain Monte Carlo (MCMC) for inference and sampled from the posterior distributions with Metropolis-Hastings and Adaptive Metropolis step methods, both available in PyMC. To speed up the model and improve the MCMC convergence, we marginalized out the latent variable N and implemented the integrated likelihood function (Guillera-Arroita et al., 2012; Royle, 2004 ). To make the integration over N values finite, we assumed 100 as the maximal possible number of individuals using a single grid cell j (25 ha). We defined the priors for the linear predictor coefficients β as diffuse normal priors N0,10-3. The priors for the measurement errors ϵi were given by N0,1/σ2 with the hyper-parameter σU0,100 (Gelman and Hill, 2006). We followed Royle et al. (2007) and we chose the gamma distributed prior G0.1,0.1 for the precision parameter of the RSR model.

The reason for not choosing a vague prior for this parameter (as e.g. in Johnson et al., 2013) was that we expected the spatial covariates included in the models would not account for all the spatial dependence alone. Part of this (unexplained) variation is likely related to species movement behavior occurring at multiple spatial scales. It is also commonly known that the variance components are poorly identified in these types of models (Royle et al., 2007). To obtain posterior distributions of parameters, we ran a MCMC sampler with three chains for 500,000 iterations each (removing the first 400.000 as a burn-in phase of the sampling process) and with the thinning parameter set to 20 to avoid autocorrelation between samples. The convergence was assessed through visual inspection of MCMC trace plots and Gelman–Rubin diagnostics provided by the PyMC software. We evaluated the fit of the models through visual inspection of standard model diagnostics plots (see Appendix 1—figures 2835). We used posterior predictive distribution and Bayesian ‘p-value’ to assess the goodness of fit of each model (Kery and Royle, 2016). The source code of our models is available at https://github.com/mripasteam/herbiscapes/ (Bubnicki et al., 2019; copy archived at https://github.com/elifesciences-publications/herbiscapes).

Ungulate community-level analysis

Request a detailed protocol

We explored the distribution of ungulates across the landscape, as predicted by species-specific models, in the context of the community. First, the spatial overlap between species was evaluated by means of pairwise Pearson correlations of their relative density surfaces. Next, we converted ungulate relative densities to herbivore biomass using the following average weights per species: red deer female 90 kg, red deer male 150 kg, roe deer 20 kg, wild boar 80 kg, moose 200 kg and European bison 400 kg (Borowik et al., 2016; Jedrzejewska et al., 1994). Total biomass was estimated as the sum of each species’ biomass for each landscape cell (25 ha pixel) in our prediction grid. We further calculated and mapped the index of the functional diversity of the entire ungulate community (FDis; Laliberté and Legendre, 2010) using the R statistical software (Core Team, 2019) and the R package FD v.1.0–12. The FDis was calculated based on the predicted relative densities of all ungulates and the following species-specific traits: body mass, diet type and gut type. The FDis is the mean distance in multidimensional trait space of individual species to the centroid of all species. It can account for species abundances by shifting the position of the centroid toward the more abundant species and weighting distances of individual species by their relative abundances (Laliberté and Legendre, 2010). The input data for a FDis calculation are 1) a matrix with species traits and 2) a matrix with relative densities that describe how much weight to assign to each individual observation. We compiled a six row (species) by three column (traits) matrix with one quantitative and two qualitative traits, namely body mass, gut type and diet type (Appendix 1—table 2).

By means of hierarchical cluster analyses (Lê et al., 2008), we grouped all landscape grid cells (25 ha each) with similar values of FDis, and total and species-specific biomass into clusters representing ecologically distinct landscape-scale herbivory regimes, or ‘herbiscapes’. Specifically, following the approach of Hempson et al. (2015), we used the R package FactoMineR v.1.39 (Lê et al., 2008) and its HCPC (hierarchical clustering on principle components) function. The HCPC requires that PCA (principal component analysis) is performed on variables prior to clustering, which limits the impact of covariance amongst variables on the subsequent clustering algorithm. To build the cluster tree, we used HCPC default values for the metric (Euclidean distance) and method (Ward’s) parameters. Similarly to Hempson et al. (2015), the number of clusters was determined by assessing the inertia (i.e. change in within cluster homogeneity) gained by cutting the tree at different levels and the ecological interpretability of the resulting clusters. Eventually, we split the cluster tree into five independent clusters (Appendix 1—figure 2).

The hierarchical clustering analysis allowed us to learn how the combined effect of predation risk and resource quality translates into the composition and abundance of the ungulate community and, in consequence, into the diversification of herbivory pressure on the ecosystem.

Vegetation analysis

Request a detailed protocol

Using an independent vegetation dataset from a large-scale inventory of tree regeneration (part of the LIFE+ ForBioSensing project, contract number LIFE13 ENV/PL/000048), we tested if the predicted variation in the landscape-scale distribution of large herbivores, synthesized into ecologically distinct herbivory regimes that is herbiscapes, affects tree browsing intensity and regenerating tree species composition. We used data collected in 2017 at 385 plots spread randomly across the entire BF. Each plot contained two concentric sub-plots: 1) with a radius of 1.3 m (area of 5 m2) at which all trees with height <30 cm excluding seedlings were recorded, and 2) with a radius of 2.52 m (area of 20 m2) at which all trees with height ≥30 cm and diameter at breast height <2 cm were recorded. Additionally, each individual tree was checked for any sign of fresh or 1-year-old browsing of its main shoot. The <30 cm tree sapling community is structured mainly by bottom-up factors (and/or forest management practices) and only minimally influenced by ungulate herbivory (Kuijper et al., 2010a), whereas the ≥30 cm tree sapling community is within the foraging height class preferred by ungulate herbivores (Kuijper et al., 2013) and is therefore largely structured by ungulate top-down factors (Kuijper et al., 2010a; Kuijper et al., 2010b).

Based on this data, for each plot we calculated the cumulative browsing intensity index, expressed as the proportion of browsed individual trees out of all tree saplings ≥ 30 cm, and the difference between the two tree height classes in the proportional shares of Carpinus betulus (Carpinus) and Acer platanoides (Acer). The latter two parameters represent a measure of recruitment from the sapling-bank (<30 cm) to the taller size-class (≥30 cm). This process is to a great extent driven by large herbivores in the studied system (see Kuijper et al., 2010a and Churski et al., 2017). We specifically focused on the response of two contrasting species, Carpinus and Acer. While both species are palatable and strongly selected by the ungulate community (see Churski et al., 2017), Carpinus is highly browsing-tolerant (a typical ‘brown-world’ species sensu Churski et al., 2017) and Acer is highly-sensitive to ungulate browsing (a typical ‘green-world’ species). Long-term exclosure studies have also shown that Carpinus typically increases while Acer decreases in dominance in response to ungulate herbivory (Kuijper et al., 2010a). As both species are very common throughout the forest, we see them as suitable indicator species for the impact of ungulate herbivory on tree species composition. We explained the variation in the calculated parameters by fitting simple linear models with two interacting factors: herbiscape and reserves. The latter was added to account for potential differences in forest structure (see Jedrzejewska et al., 1994) and ungulate behaviour (see Kamler et al., 2008) between protected and unprotected areas in BF.

Results

Large carnivores

The main factor associated with the space use of both large carnivore species was human activity, as indicated by the positive effect of distance to major settlements on their spatial distributions (Figure 3). For the wolf, the density of protected areas was another important variable related to its landscape use. Wolves more often used large nature reserves (including BNP) than parts of the landscape dominated by managed forest. There was no statistically important effect of elevation, distance to touristic trails and density of touristic infrastructure on landscape distribution of neither the wolf nor lynx. However, wolves tended to use lower areas more intensively. For the lynx there was a clear tendency to use parts of the landscape further away from major human settlements; however, this effect was less evident than for wolf (the credible intervals overlapped at zero, Figure 3D). The density of protected areas had no effect on lynx distribution.

Spatial distribution of large carnivores is determined by human activity.

(a, d) Estimated parameters for the wolf and lynx models, (b, e) maps of the predicted variation in their landscape use (relative density surfaces) and (c, f) fitted spatial random effects (SRE). Specifically, panels (b) and (e) present the spatial predictions for the parameter λ, which is the expected number of individuals using a given landscape grid cell (25 ha pixel) during the sampling period (see the general and formal description of the model in the Materials and methods section). The fitted values of SRE are deviations from 0 at log-scale. The SRE captured all spatial variation not explained by the covariates included in the model, indicating parts of the landscape for which higher (“hot-spots”) or lower (“cold-spots”) activities of species were observed in the data than was predicted by the “fixed” part of the model. On panels (a) and (b) we additionally marked in bold the credibility intervals (quantile based) at which the estimated parameter values differed from 0 (for the level 2.5–97.5%, there is a 0.95 probability that the true parameter value lies within this range). The full posterior distribution for each parameter, together with its numerical description can be found in Appendix 1—figures 2627).

https://doi.org/10.7554/eLife.44937.005

To test the quality of our predictions, we compared them with existing radio-tracking data from collared wolves collected over 20 years ago (1994–1999) in BF. This showed that our model fitted to the camera trapping dataset (471 wolf detections) not only conforms to the general pattern of the wolves’ space use determined by telemetry 20 years ago, but also reflects the wolves’ response to ongoing environmental changes in the study area. See the Appendix 1 and Figure 4 for more details.

Ungulates – landscape use

In contrast to the other ungulates, the red deer, the main prey of the wolf (Jędrzejewski et al., 2002), was the only species whose landscape use was associated with that of wolves (Figures 4 and 5). The predicted relative densities of both red deer females and males were lowest in parts of the landscape intensively used by wolves, and this effect was more pronounced for females. Lynx distribution was not related to the landscape use of any ungulate species. Red deer females were also positively associated with protected areas (BNP and nature reserves) and mixed deciduous forest with an intermediate level of landscape openness (Figure 4). Red deer males showed similar tendencies, but these factors were not statistically important predictors of their landscape distribution. Instead, red deer males showed a negative association with distance to human settlements, which could have resulted from their more intensive use of open meadows surrounding settlements, especially during the rutting period. Forage habitat availability was the main factor associated with wild boar distribution, as indicated by its clear association with closed-canopy and deciduous species dominated forest stands (as indicated by the strong negative effect of coniferous species dominated forest stands; Figure 4). The presence of deciduous forest was also the main factor positively associated with the distribution of bison, followed by proximity to human settlements (as indicated by the negative association with distance to human settlements). The latter is likely related to the presence of meadows surrounding settlements, which provide optimal foraging habitats for bison (Bocherens et al., 2015; Cromsigt et al., 2012). The moose was the only species whose landscape-scale distribution was positively associated with low elevation areas (river valleys and wetlands). Roe deer were positively associated with intermediate levels of landscape openness. However, no other predictors that we used were related to the distribution of roe deer, which may be due to the low density of this species in our study system. All ungulate species except red deer males showed some level of remaining spatial auto-correlation not explained by the raster covariates included in the models (Figure 6). The spatial random effects were particularly large for the bison, whose spatial distribution in BF is strongly influenced by supplementary winter feeding and use of open areas outside the forest throughout the year (Kowalczyk et al., 2011). An interesting hot-spot of unexplained variation in distribution of bison and red deer females was the south of the national park, in an area without hunting and rich old-growth deciduous stands and with known higher densities of deer (Jędrzejewska et al., 1997).

Spatial distribution of each ungulate species was related to a unique combination of bottom-up and top-down factors.

(a) Estimated effects of covariates explaining the spatial variation in the parameter λ, which is the expected number of individuals using a given landscape grid cell (25 ha pixel) during the sampling period (i.e. the relative density; see the general and formal description of the model in the Materials and methods section). In bold, we marked the credibility intervals (quantile based) at which the estimated parameter values differed from 0 (for the level 2.5–97.5%, there is 0.95 probability that the true parameter value lies within this range). The full posterior distribution for each parameter, together with its numerical description can be found in Appendix 1. (b) For easier interpretation of the effects of landscape openness and percentage share of coniferous species (together with their quadratic terms and assuming an average level for the other covariates), we plotted predictions of the relative density for each species and for each variable. The values of the covariates were scaled and zero-centered before making predictions; thus 0 corresponds to the average value of a given covariate. The full posterior distribution for each parameter, together with its numerical description can be found in Appendix 1—figures 1425).

https://doi.org/10.7554/eLife.44937.006
Only red deer (the major wolf prey) showed a negative spatial association with parts of the landscape intensively used by wolves, whereas distributions of other species were shaped by bottom-up factors.

The maps show substantial variation in predicted landscape use for the five studied ungulate species (relative density surfaces). Specifically, this figure presents the spatial predictions for the parameter which is the expected number of individuals using a given landscape grid cell (25 ha pixel) during the sampling period (see the general and formal description of the model in the Materials and methods section). The predictions are based on the parameter estimates presented in Figure 3 and include the spatial random effects (see Figure 4).

https://doi.org/10.7554/eLife.44937.007
All ungulate species except red deer males showed some level of remaining spatial auto-correlation not explained by the raster covariates included in the models.

The maps show the fitted spatial random effects (SRE) for the models of five studied ungulate species. The SRE are deviations from 0 at log-scale. The SRE captured all spatial variation not explained by covariates included in the model, indicating parts of the landscape for which higher (‘hot-spots’) or lower (‘cold-spots’) activity of species was observed in the data than predicted by the ‘fixed’ part of the model. The SRE were particularly large for the bison, whose spatial distribution in BF is strongly influenced by supplementary winter feeding and the use of open areas outside the forest during the year.

https://doi.org/10.7554/eLife.44937.008

Ungulates – detection rates

At the scale of single camera trap sites (1 ha), increased canopy openness led to higher detection rates of red deer and moose (Figure 7). This could be explained by cameras having better detection in open areas (Marcus Rowcliffe et al., 2011) and/or a preference of these species to forage in canopy gaps (Churski et al., 2017; Kuijper et al., 2009). Red deer females were associated with even larger canopy gaps, whereas for red deer males and moose there was an optimum value of canopy openness, above which the detection rate started to decrease. Roe deer followed a similar pattern as red deer but there was no statistical support for this result. At-site detection rate of wild boar was highest in deciduous forest patches with relatively closed canopies (in line with Kuijper et al., 2009). However, wild boar, as well as red deer, roe deer and bison were also relatively frequently detected at forest patches dominated by coniferous tree species, indicating a context-dependence in the selection of small habitat patches. For example, as forage availability changes between seasons, deciduous patches may be preferred in the green season while coniferous patches in winter. Temperature affected detection rate non-linearly by influencing the activity level of red deer males, with the optimum at the yearly average temperature. In the case of the wild boar the highest detection rates were found at high temperatures, coinciding with their reproductive period in summer. And in case of the moose, the highest activity level was found at moderately high temperatures. The bison was the only species whose detection rate was strongly affected by snow cover. Winter supplementary feeding causes bison to aggregate near feeding stations or outside BF (Kowalczyk et al., 2011), hence dramatically decreasing their detection rate in other parts of the forest.

Factors influencing detection rates of ungulate species.

(a) Estimated effects of covariates explaining the variation in the parameter γ, which is the expected (daily) detection (trapping) rate at a single camera trap site (see the general and formal description of the model in the main text in the Materials and methods section). The spatial covariates were calculated at a pixel resolution of 100 m (1 ha) and the temporal covariates were calculated as the average value for the whole period of camera trapping at a given location. In bold, we have marked the credibility intervals (quantile based) at which the estimated parameter values differ from 0 (for the level of 2.5–97.5 %, there is a 0.95 probability that the true parameter value lies within this range). The full posterior distribution for each parameter, together with its numerical description can be found in Appendix 1. (b) For easier interpretation of the effects of canopy openness, the percentage share of coniferous species and temperature (together with their quadratic terms and assuming an average level for other covariates), we plotted predictions of the detection rate for each species and for each variable. The values of the covariates were scaled and zero-centered before making predictions; thus 0 corresponds to the average value of a given covariate. The full posterior distribution for each parameter, together with its numerical description can be found in (Appendix 1—figures 1425).

https://doi.org/10.7554/eLife.44937.009

Ungulates – spatial overlap

All ungulate species except moose showed some level of pairwise positive spatial associations as indicated by the Pearson correlation of their relative density surfaces predicted by the models (Figure 8). However, the strength of these associations was relatively low, indicating substantial spatial variation in the structure of the whole community. Unsurprisingly, the strongest overlap in space was between red deer females and males. However, the estimated correlation (0.71) was far from a ‘perfect’ overlap, indicating red deer are sexually segregated in space in BF (see Kamler et al., 2008). The relatively high spatial overlap between red deer males and roe deer (0.58) likely resulted from their utilization of similar parts of the landscape, often close to the forest edge and large clearings with human settlements. The moose was a clear exception showing a negative spatial association with all other ungulate species. This indicates that the moose has a specific (spatial) niche in our system, strongly associated with low lying areas like river valleys and wetlands (see Figures 2, 4 and 5). Interestingly, when comparing the surfaces of spatial random effects, there was a relatively strong positive spatial association (0.59) between bison and wild boar (Figure 6). A possible explanation for this pattern may be that wild boar are attracted to supplementary food at bison feeding stations, where next to hay and silage, beetroots are provided for bison (Kowalczyk et al., 2011).

There was substantial spatial variation in the structure of the ungulate community, as indicated by the weak spatial associations between species.

The pairwise spatial overlap between the studied ungulate species is presented as the Pearson correlation heatmap of their relative density surfaces predicted by the models (a) and fitted spatial random effects (b).

https://doi.org/10.7554/eLife.44937.010

Functional diversity and hierarchical clustering analysis

When mapped, both the estimated total biomass of the ungulate community and the functional diversity index (FDis) showed distinct patterns across the studied landscape (Figure 9). Interestingly, the lowest values of both parameters were associated with coniferous-dominated tree stands at higher elevations (Figures 2 and 5) belonging to the least productive parts of this landscape (Faliński and Falińska, 1986; Kwiatkowski, 1994). By means of hierarchical cluster analyses, we further grouped all landscape grid cells (25 ha each) with similar values of FDis, and the total and species specific biomasses. We identified five clusters, characterized by different sets of risk-related and environmental factors and composed of different sets of ungulate species (Figure 10). The ‘red’ cluster (id = 1) was characterized by high wolf and lynx use and low quality foraging habitat (low elevation, mixed tree stands with large shares of coniferous tree species, Figure 10B). In terms of ungulate biomass, this cluster was dominated by moose, red deer and wild boar (Figure 10C). However, red deer relative density was the lowest out of all five clusters. The total ungulate biomass in this cluster was low but the FDis value was relatively high (Figure 10D). A large part of the ‘red’ cluster was spatially associated with the lowest elevated areas, that is marshlands and river valleys (compare Figure 10A with Figure 2). The ‘green’ cluster (id = 3) was characterized by high predator presence and high-quality foraging habitat (moderate elevation, mixed tree stands with large shares of deciduous tree species), which increased the numbers of seemingly risk-insensitive species like wild boar and bison. Red deer males, which had a less pronounced negative spatial association with wolves than females (Figure 4), also had a higher biomass in this cluster than in the ‘red’ one. Moreover, the ‘green’ cluster was characterized by low values of total biomass and high values of FDis. This cluster in the landscape, mainly covered remote areas of low human activity and high use by wolf, and occurred in protected areas (BNP and nature reserves). The ‘blue’ cluster (id = 2) was characterized by moderate use by predators and low-quality foraging habitat (high elevation, coniferous dominated tree stands), and had the lowest values of both total biomass and FDis. However, the wolf used these areas less intensively than the ‘red’ and ‘green’ clusters, and the biomass of red deer was higher than in those high wolf use clusters. The ‘purple’ cluster (id = 4) was characterized by low use by predators and high-quality foraging habitat (high elevation, mixed tree stands with large share of deciduous tree species), and was dominated by red deer. Both sexes of red deer were at their most abundant in this cluster, which covered areas of high human activity and the lowest use by wolf. These seemingly ‘safe’ parts of the landscape were also characterised by high-quality habitat for red deer – a mosaic of deciduous and mixed tree stands at higher elevations and relatively close to forest edges and large clearings surrounding human settlements. The roe deer was also spatially associated with this cluster. The ‘orange’ cluster (id = 5) was characterized by moderate use by predators and high-quality foraging habitat (moderate elevation, deciduous dominated tree stands) and had the highest values of total biomass and FDis. This cluster was dominated by the bison - the main grazer and largest species in our system. Also, the wild boar was at its highest relative density in this cluster. It is worth mentioning that a large part of this cluster was in the southern part of Białowieża National Park, comprising some of the best preserved parts of this forest.

Maps of estimated total biomass [kg/pixel] and functional diversity (FDis).

Both variables, together with species-specific biomass raster layers, were used in the hierarchical clustering analysis. Interestingly, the lowest values of both parameters were largely associated with coniferous-dominated tree stands at higher elevations (compare with Figure 2), often re-planted after extensive clear-cuts from the past and belonging to the least productive parts of this landscape (Faliński and Falińska, 1986; Kwiatkowski, 1994).

https://doi.org/10.7554/eLife.44937.011
Spatial variation and functional diversity of a herbivore community within a temperate forest herbivome, decomposed into landscape-scale herbivory regimes – herbiscapes.

(a) The spatial distribution of the five identified clusters (‘herbiscapes’) based on hierarchical clustering on the principal components analysis. (b) The distribution of identified clusters in the space of the first two principal components. (c) Tukey-like boxplot comparing the species-specific and total biomass in each of the identified clusters. (d) Tukey-like boxplot comparing the functional diversity index (FDis) of the ungulate community between the identified clusters. The boxplots show median values and the lower and upper hinges correspond to the first and third quartiles. The colors of clusters are consistent amongst all figures. Each data point on panel (b) corresponds to a 25 ha pixel mapped on figure (a). On panel (b), for ease of interpretation, we additionally plotted (blue arrows) the major environmental and risk-related gradients in the study area which were used as covariates in the ungulate models.

https://doi.org/10.7554/eLife.44937.012

Spatial patterns in herbivores and impact on woody vegetation

The highest proportions of browsed trees were found in herbiscapes 3, 4 and 5 (Figure 11c), which are characterized by high values of total ungulate biomass and/or high functional diversity index (FDis). These three herbiscapes covered the more fertile parts of the landscape dominated by productive deciduous and mixed forests (Figure 10B) with abundant tree regeneration dominated by Acer platanoides and Carpinus betulus (Acer and Carpinus hereinafter, Figure 11d). The variation in browsing intensity between the herbiscapes was reflected in the recruitment patterns towards taller tree sapling size classes (>30 cm) of Acer and Carpinus. In accordance with previous experimental exclosure studies in BF (Churski et al., 2017; Kuijper et al., 2010a; Hedwall et al., 2018), the proportion of the palatable but browsing intolerant Acer in the community of regenerating trees decreased as herbivore pressure increased (Figure 11a,d). In contrast, the palatable but highly browsing-tolerant Carpinus, showed the opposite pattern (Figure 11b,d). These supposedly herbivore-driven shifts in the species composition of regenerating trees were most pronounced in nature reserves, whereas areas outside the reserves broadly showed the same patterns but less closely followed the observed patterns in cumulative browsing intensity.

The patterns of browsing intensity and the composition of regenerating tree species differ between herbiscapes and between the protected and managed parts of the BF landscape.

(a, b) The mean differences between the two tree height classes (saplings < 30 cm and ≥30 cm) in the proportional shares of Acer platanoides (Acer) and Carpinus betulus (Carpinus). The numbers on top of both figures are the proportions of plots without a single individual of each species. (c) The cumulative browsing intensity index, expressed as the proportion of browsed individual trees out of all tree saplings ≥ 30 cm. The numbers on top of this figure are the proportions of plots with no browsing. Panels (a), (b) and (c) show the values predicted by the fitted linear models with two interacting factors: herbiscape (six levels) and reserves (two levels coded as 0: unprotected areas and 1: protected areas). The error bars are standard errors of the mean. (d) The mean proportions of selected tree species (saplings < 30 cm and ≥30 cm) out of the entire pool of the main regenerating trees recorded at the sampled plots, showed for the herbiscapes and protected (R1) and unprotected (R0) areas separately. The value on top of each stacked bar is the number of plots used to calculate the mean proportion. Cb: Carpinus betulus, Ap: Acer platanoides, Tc: Tilia cordata, Qr: Quercus robur, Pa: Picea abies, Ot: other species (Betula pendula, Sorbus aucuparia, Fraxinus excelsior, Alnus glutinosa, Populus tremula, Pinus sylvestris, Ulmus).

https://doi.org/10.7554/eLife.44937.013

Discussion

Our results show how variation in spatial distribution in a community of large herbivores in a temperate forest ecosystem can be explained by species-specific associations with major ecological gradients operating at the landscape scale. This creates ecologically distinct landscape-scale herbivory regimes (‘herbiscapes’), which are interactively driven by both environmental bottom-up and biotic top-down (large carnivores) factors in combination with human-driven (cascading) effects. In addition, our analyses suggest that these herbiscapes differ in browsing intensity and impact on vegetation, indicated by the changing proportions of recruitment of browsing-sensitive versus browsing-tolerant tree species (Figure 12).

Graphical abstract presenting the steps taken to reveal the relationships uncovered in this study.

Image credit: Lisa Sánchez Aguilar

https://doi.org/10.7554/eLife.44937.014

Spatial heterogeneity in the landscape distribution of large herbivores

All ungulates in our study showed specific, non-uniform distributional patterns that were associated with species-specific combinations of bottom-up and/or top-down forces, including predation, human presence, availability of resources and the (bio)physical properties of the landscape (Figures 4 and 5). The substantial spatial variation in each of these landscape components, when combined, resulted in an aggregated, non-uniform pattern of landscape distribution of all ungulates, as predicted by Fryxell (1991) and Hopcraft et al. (2010). The spatial distribution of each species (Figure 5) followed a characteristic shape of a hollow or sigmoidal curve (Figure 13) when plotted as a graph of ranked relative densities predicted for each 25 ha landscape pixel. This pattern of spatial variation in abundance of different populations has been shown to be universal over large spatial domains and for different taxa (Brown et al., 1995). Our study contributes to this finding by showing that similar spatial patterns can be observed within a population and within a local landscape. The shape of these ranked-abundance curves is informative for the properties of a given continuous (landscape) surface (Rocchini and Neteler, 2012), for example the sigmoidal curve of wild boar indicates a highly heterogeneous distribution with many hot and cold spots (high and low use), whereas the hollow curve of moose indicates a rather homogeneous distribution at low density with only some hot-spots.

Ranked relative densities predicted for all ungulate species and for each 25 ha landscape pixel.

All plots follow a characteristic shape of a hollow or sigmoidal curve (Brown et al., 1995; Rocchini and Neteler, 2012) indicating a non-uniform, heterogenous distribution of all studied species.

https://doi.org/10.7554/eLife.44937.015

Interestingly, the degree of spatial aggregation of different ungulate species in our study system was not related to their body size. This contrasts with studies from African savannas where larger herbivore species were more evenly distributed over the landscape than the smaller species (Cromsigt et al., 2009; du and Owen-Smith, 1989). The two largest herbivores in our study system, the European bison and moose, showed highly aggregated landscape distributions. The moose was strongly associated with lower elevation habitats (i.e. river valleys and wetlands), whereas landscape use by bison was strongly related to locations of supplementary winter feeding and open areas inside the forest, which are distributed sparsely within the study area and often associated with anthropogenic activity. The distribution pattern of bison could also be the result of partial migrations outside the forest in winter (Kowalczyk et al., 2013; Kowalczyk et al., 2011). This mainly human-driven distribution of bison was also manifested in the largest values for unexplained spatial variation in the model output (fitted spatial random effects, Figure 6).

The distributions of all species, when combined, revealed substantial spatial variation in the composition of the entire ungulate community (Figure 10). To our knowledge this is the first empirical study presenting a synthesized, high-resolution and spatially explicit approach (sensu Royle et al., 2007) combining bottom-up and top-down factors to explain the landscape-scale variation of an entire large herbivore community in a temperate forest ecosystem. Including both resource- and predator-related factors was critical for achieving this goal as they can operate simultaneously and interactively (Anderson et al., 2010; Fryxell, 1991; Hopcraft et al., 2010). Similar studies in African ecosystems have shown that the distribution and diversity of a community of large herbivores can be driven by bottom-up (e.g. habitat heterogeneity; Cromsigt et al., 2009), top-down (e.g. anthropogenic fire; Klop and Prins, 2008) or interactive effects of both factors (e.g. distance to water and settlements; Ogutu et al., 2010). More recent studies, also from African ecosystems, have extended this approach by including predation-related factors, revealing the trade-offs native ungulates make to cope with changes in forage availability, human disturbance and predation risk (Schuette et al., 2016) and showing that a top predator can have species-specific spatial associations with herbivores (Anderson et al., 2010). For example, the latter study showed that lions were positively associated with large-bodied migratory ungulates but negatively associated with smaller non-migrants. Most comparably to our study, Anderson et al. (2016) used an extensive network of camera traps and a spatially-explicit occupancy modeling framework to quantify the spatial distribution of African savanna herbivores. Interestingly, they quantified pairwise interactions between all modeled species demonstrating the emergence of strong positive spatial associations among a diverse group of savannah herbivores. This is in line with our results where all ungulates except moose showed positive spatial associations. However, in our study, these spatial associations were rather weak, indicating substantial spatial variation in the structure of the ungulate community. Our contribution to all the above studies is a high-resolution picture of the spatial structure of an entire community of large herbivores that incorporates both bottom-up and predation- and human-related top-down factors. We believe one of our major points of novelty is in providing information on these spatial interactions for a temperate ecosystem.

The effect of large carnivores on the spatial structure of the large herbivore community

The red deer was the only species whose landscape use was associated with that of large carnivores, with lower red deer presence in areas with higher wolf use. This is in line with other studies that have shown that predator top-down effects can work selectively on some members of large herbivore communities (Sinclair et al., 2003; Valeix et al., 2009) and have a complex and context-specific nature especially in multiple-prey and multiple-predator systems (Davies et al., 2016; Moll et al., 2016). However, wolves, by seemingly affecting the spatial distribution of red deer (see below), one of the dominant species in BF, re-structured the composition of the entire community of large herbivores (Figure 14) and increased the degree of its spatial heterogeneity (Figure 10).

Scatter plot showing how wolves, by affecting the spatial distribution of red deer, re-structured the composition of the entire community of large herbivores.

The predicted biomass of red deer (both sexes, X axis) is shown as a share of the total biomass of the entire community of ungulates. Each dot on the plot is a 25 ha landscape pixel. In color, we marked relative wolf encounter risk (scaled values of the parameter λ) as predicted by the model fitted to the wolf camera trapping data. We interpreted the parameter λ as large carnivore space use intensity and assumed that from a prey perspective, λ is proportional to the predator encounter rate. The dashed line is a reference line indicating the y=x relationship.

https://doi.org/10.7554/eLife.44937.016

The red deer is the main prey species for wolves in our study area and experiences the highest predation pressure by wolves, in contrast to the European bison, moose and roe deer, which comprise only a small proportion of their diet, and wild boar, which is a secondary prey species (Jędrzejewska et al., 1997; Jędrzejewski et al., 2002). Hence, this might explain why the the red deer was also the species that seemed to react most strongly to the space use of its main predator. The most intensively used part of a wolf territory (on an annual basis) is related to the location of their dens during the reproductive period, and is generally far from human settlements (Jędrzejewski et al., 2001; Kuijper et al., 2015). These high wolf use areas likely only have higher predation rates during the reproductive period (Kuijper et al., 2013). Over the rest of the year, wolves move more widely across their territories (Jędrzejewski et al., 2001) and kills may be distributed more widely across annual wolf territories. Hence, the areas that seemed to be avoided by red deer are not necessarily areas with the highest predation risk on an annual basis. This finding adds to the growing recognition that prey species perceive risk based on various factors such as the space use of large carnivores and physical landscape and not necessarily by kill site distribution (Kohl et al., 2018; Gaynor et al., 2019).

The distributional pattern of red deer was similar for both sexes (one independent model was fitted for each), but females had a stronger negative association with wolf space use than males. This indicates that females are more sensitive to wolf presence than males and is consistent with the selective killing of this sex and juveniles by wolves in our study system (Jędrzejewski et al., 2000; Jędrzejewski et al., 2002). These apparent effects of wolf space use on red deer distribution could have resulted from both non-lethal (behaviorally mediated) as well as lethal (density-mediated) effects. With our data we were not able to distinguish between these two mechanisms, although previous studies have shown that non-lethal risk effects play an important role in affecting the responses of red deer at both fine- and large-scales in this system (Kuijper et al., 2015; Kuijper et al., 2013; Kuijper et al., 2014).

In contrast to the red deer, environmental bottom-up factors, particularly landscape topography and resource availability (natural or supplemented by humans), had the strongest associations with the spatial distribution of the other, less predation-sensitive ungulate species (Figures 4 and 5). Although the wild boar is a secondary prey species for wolves (Jędrzejewski et al., 2002), its high abundance during t study period means that predation by wolves contributed little to its annual mortality (Jędrzejewska et al., 1997). This may explain why we did not observe a negative association between wild boar and wolf space use. Moreover, a previous study of ours found that the wild boar displayed behaviour suggesting it perceived no predation risk in response to the presence of fresh wolf or lynx scats (Kuijper et al., 2014; Wikenros et al., 2015).

Our results seem contrary to those of Theuerkauf and Rouys (2008), who carried out a similarly focused study in the same study area with use of pellet counts. They concluded that habitat alteration by forest exploitation and hunting by humans influenced the density distribution of ungulates, including red deer, more than predation risk by wolves. Despite our results indicating a spatial mismatch between the red deer and wolves’ landscape use, this is likely to be caused by a cascading effect of humans. Moreover, the patterns of ungulate space use shown by our analyses revealed an additional source of variation, which involves a species-specific response to the inter-relationships between human and predator space use.

We were surprised to find that the Eurasian lynx – the other large carnivore in our study area – had no apparent effect. As ungulates constitute the bulk of the lynx’s diet, with the roe deer being the major prey (60% of the diet) and the red deer being the alternative prey (22%) (Okarma et al., 1997), it is an important apex predator in this system. Lynx predation is a major mortality factor for both cervids, taking 21–36% of roe deer and 6–13% of red deer population numbers (Okarma et al., 1997; Jędrzejewski et al., 1993; Jędrzejewska et al., 1997; Jędrzejewska and Jędrzejewski, 2005). Moreover, it was recently revealed that red deer clearly react with anti-predatory behavior to olfactory cues of the lynx in BF (Wikenros et al., 2015). It is thus striking that these highly sensitive behavioral responses do not lead to changes in the spaces use of ungulates at the landscape scale. This may be a result of the combined effects of the low densities of both the roe deer and the lynx and the particular hunting mode of lynx, which, typically for felids, relies on fine-scale habitat characteristics that allow the predator to exploit prey independently of their spatial distribution (Podgórski et al., 2008; Schmidt, 2008). Although there are many studies on the impact of large carnivores on the space use of their prey species in temperate systems, the majority of these have focused on single carnivore - single prey relationships (see e.g. Creel et al., 2005; Kauffman et al., 2007; Lima and Dill, 1990; Mao et al., 2005). Our study is the first we are aware of to show how these carnivore top-down effects structure an entire ungulate community in a temperate landscape. This knowledge is relevant as most terrestrial ecosystems are characterized by a diversity of herbivorous prey species and the differential responses of functionally similar prey species to apex predators can to a large extent determine the potential for trophic cascading effects of large carnivores (Ford and Goheen, 2015; Ford et al., 2015; Rosenfeld, 2002).

Mediating predator-prey interactions at a landscape scale by humans

Humans can drive complex interactions between species, particularly by affecting keystone species like large carnivores (Worm and Paine, 2016). We found that humans were the main factor associated with the spatial distributions of both the wolf and lynx, which had lower activities in parts of the landscape heavily used by humans (in line with Theuerkauf et al., 2003). In this way, human presence can be beneficial for ungulate prey species as large carnivores generally avoid human presence and activity more strongly than their ungulate prey species (Rogala et al., 2011), leading to so-called ‘human shields’ (Berger, 2007). The observed spatial (re)distribution of red deer in our system seems to result from human-induced shifts in space use of the wolf. These kinds of three-way trophic cascades (following the definition of Ripple et al., 2016) involving humans, large carnivores and ungulates have been found in different systems, mainly in landscapes with moderate human activity (Berger, 2007; Hebblewhite et al., 2005) and are likely much more pronounced in highly human-dominated landscapes (Kuijper et al., 2016). We carried out our study in the forest considered to be the best preserved lowland forest ecosystem in Europe, which is replete with natural wildlife communities and ecological processes (including trophic interactions) still operating at the landscape scale. Thus, we believe it provides new knowledge on the fundamental structuring and functioning of European forest ecosystems that will help to predict the ecological effects of the ongoing large carnivore recolonisation of more human-altered habitats across Europe (Chapron et al., 2014).

Landscape scale herbivory regimes and their potential consequences for vegetation

How the consumptive off-take by large mammal herbivores varies in space at a continental-scale has recently been shown by Hempson et al. (2015), who classified African herbivore communities into ecologically distinct herbivory regimes, or herbivomes. In our study, we down-scaled this approach to the landscape level, exploring the spatial variation and functional diversity of the herbivore community within a single herbivome, decomposing it into ecologically distinct landscape-scale herbivory regimes. We refer to these as herbiscapes. Additionally, we extended the herbivome concept by including the direct and indirect top-down effects of higher trophic levels: large carnivores (keystone species) and humans (hyper-keystone species). Thus, similarly to a spatial ‘profile’ of each single species, the major dimensions of a herbiscape are oriented along both environmental (landscape feature, resource availability) and trophic interaction induced (direct and perceived predation risk, human activity) gradients operating over a landscape (Figure 10).

As large herbivores influence many important processes of ecosystem functioning (Hobbs, 1996), the observed landscape-scale variation in the structure of the large herbivore community (i.e. herbiscapes), when stable enough, could result in a differential impact of herbivores on (woody) vegetation. Our tree regeneration analysis revealed that the patterns of ungulate browsing intensity followed the variation in the modeled community-level distribution of ungulates in BF. The highest browsing intensities occurred in herbiscapes characterized by high ungulate biomass and/or functional diversity. Compositional changes in regenerating trees across herbiscapes were in accordance with our earlier experimental studies (Churski et al., 2017; Kuijper et al., 2010a): with higher levels of browsing intensity the proportion of a common browsing-tolerant tree increased and the proportion of a browsing-intolerant tree decreased. These patterns strongly suggest that the impact of ungulate communities on tree regeneration varies greatly at the landscape-scale, even in relatively homogeneous forest landscapes like BF. Earlier experimental studies showed that the species composition of the small tree sapling community in BF is shaped by environmental bottom-up factors, whereas herbivores are the main factor limiting recruitment towards taller size classes (above 50 cm; Kuijper et al., 2010a). We therefore interpret the observed shift in dominance towards browsing-tolerant species in the taller height class (>30 cm) as the result of herbivore top-down effects, as indicated by the clear differences in browsing intensity (Kuijper et al., 2010a; Churski et al., 2017). Humans are an inherent part of these multi-trophic interactions as they affect both carnivore and species-specific herbivore space use, cascading down in a complex but predictable way to the lower trophic levels. Our study suggests, that the existence of more stable herbiscapes at larger spatial scales could create a mosaic of differential impact on woody plant communities and consequently spatial variation in expressed tolerance traits. Herbiscapes could therefore strongly contribute to the creation of a heterogeneous patchwork of green (less browsing tolerant) and brown (more browsing tolerant) worlds (sensu Churski et al., 2017) driven by variation in herbivory pressure at multiple spatial scales in temperate forest systems.

Recent views on the biome concept have emphasized the role of biotic interactions (e.g. herbivory) in creating multiple biome stable states under similar climatic conditions (Moncrieff et al., 2016; Woodward et al., 2004). This approach has proven to be very useful in explaining the spatial distributions of different vegetation communities maintained by functionally distinct guilds of herbivores at continental-scales (Charles-Dominique et al., 2016; Hempson et al., 2015). On the basis of the high-resolution data we collected on ungulate species distribution, our study revealed the presence of a large continuous variation in herbivore community structure at a much finer, within-landscape scale. Although European assemblages of wild herbivores are not as species-rich as those found in African savannas, the species we studied are clearly functionally diversified (see for example Hofmann, 1989). These differences in resource use and foraging behavior among species within the community could differentially impact the vegetation if certain parts of the landscape are dominated by different functional types of herbivores for a long enough period.

Studies from African savannas have indicated that the stability of these relationships is poorly understood, but is likely maintained by strong large-scale positive feedbacks between vegetation, abiotic resources and consumers (Charles-Dominique et al., 2016). This long-term stability likely provides the basis for co-evolutionary dynamics between functional types of herbivores and plants in these systems. Whether such feedbacks could also be operating at smaller spatial scales (i.e. in herbiscapes within a herbivome) and play a role in creating small-scale variation in vegetation structure is an intriguing question. While we do not believe they could be stable on an evolutionary time-scale, there could be sufficient stability over several decades to significantly impact woody plant communities (Kuijper et al., 2010a). In our system, these feedbacks have been controlled by landscape scale anthropogenic factors that have been present in a similar spatial arrangement for decades or even centuries (i.e. villages and roads have their origins in historical times; see for example Samojlik et al., 2016). As a result, the distribution of herbiscapes, driven by human-induced carnivore space use is arguably stable enough to create differential impact on vegetation in different herbiscapes. As humans are a crucial factor determining and restricting space use of large carnivores worldwide (Tucker et al., 2018), they likely contribute toward stabilizing herbiscapes in many human influenced landscapes across the globe.

Similarly, the results of our study have implications for patterns in vegetation structure that have been observed at large spatial scales. Herbiscapes will likewise also occur within the herbivomes described for the African continent. Our study indicates that predators can be a main factor structuring herbivore communities by redistributing predation-sensitive species over the landscape. As African systems harbor one of the most carnivore-rich communities in the world (Ripple et al., 2014), large carnivores are expected to greatly influence variation in herbivore community structure within African landscapes (Thaker et al., 2011; Valeix et al., 2009). A major difference between our, and African study systems, is that humans supposedly play a less pronounced role in determining the space use of large carnivores and herbivores (and hence the spatial arrangement of herbiscapes) in the latter (but see Tucker et al., 2018). The importance of this landscape-scale (or within-herbivome) variation in herbivore community structure in creating clear spatial patterns in vegetation, have already been illustrated for African systems. For example, Ford et al. (2014) showed that both predation risk and plant defenses enabled plants to thrive in different parts of a landscape; consequently, the thorniness of tree communities decreased across the landscape, contributing to intra-biome variation driven by predator-prey interactions (Ford et al., 2014).

The knowledge obtained in this study on spatial variation in the densities of local wildlife populations dependent on species-specific responses to habitat and disturbance factors gives us a better understanding of wildlife communities and may be relevant for effective wildlife management. In contrast to the common assumption that for management purposes wildlife populations are uniformly distributed, our study emphasizes the existence of large spatial variation in the landscape scale densities of ungulates. Such an approach fits particularly well in the recently developed concept of ‘hunting for fear’, which promotes the spatio-temporal diversification of management techniques based on perception of risk in wildlife (see Cromsigt et al., 2012). Moreover, maps visualizing cold- and hot-spots of ungulates across landscapes, together with information about community compositions, could also be useful tools for wildlife managers.

Our results also introduces valuable knowledge relevant to the conservation of natural habitats. We argue that a unique, and unfortunately still undervalued character of Białowieża Forest is that ecosystem processes are still operating here at a landscape scale, as shown by our study. The ‘herbiscapes’ proposed in this paper are another unique aspect of the ecological processes that should be preserved within this area and will contribute to the ongoing debate on future conservation strategies for Białowieża Forest as UNESCO World Heritage (Mikusiński et al., 2018).

Conclusions

In conclusion, our study has illustrated that spatial variation in the structure of an entire large herbivore community results from interactive effects of species-specific responses to major ecological gradients operating at the landscape scale. Humans were a crucial factor associated with the landscape use of wolves and lynx. While lynx were not associated with the space use patterns of any ungulate, wolves were strongly (negatively) associated with the spatial distribution of their main prey species (red deer), affecting the ungulate co-occurrence patterns at the landscape scale. The space use of European bison, moose, roe deer and wild boar was related to food resources. These processes led to the distribution of different functional types of herbivores over the landscape and created a clear spatial structure in the herbivore community, which we referred to as herbiscapes. Vegetation analyses suggested that herbivore impact measured by browsing intensity and regeneration of browsing-tolerant tree species consistently differed between herbiscapes. When these herbiscapes are stable enough they could be an important mechanism driving variation in herbivore impact on woody vegetation and thus maintain heterogeneity in a wide range of ecosystems.

Appendix 1

Linking spatial patterns of terrestrial herbivore community structure to trophic interactions

Jakub W. Bubnicki, Marcin Churski, Krzysztof Schmidt, Tom A. Diserens, Dries J.P. Kuijper

Methods

Spatial covariates – data sources and processing

The GIS layers used as spatial covariates in our models were prepared using QGIS (QGIS Development Team, 2017) and GRASS GIS (Neteler et al., 2012) open source software. We compiled a set of high-resolution (5 m) rasters to describe environmental and human-related gradients within the study area. All raster covariates were summarized at two focal resolutions: 100 m (camera trap site scale) and 500 m (landscape scale). The density-based covariates, that is density of protected areas and density of touristic infrastructure were calculated for each pixel using 5 × 5 km kernel windows. The density of protected areas included Białowieża National Park and all strictly protected nature reserves in the study area. The density of tourist infrastructure was calculated based on locations of so-called ‘Points of Interest’ (POIs; e.g. hotels, restaurants, agrotourism farms, museums, visitors centers, camp sites but also smaller touristic infrastructure like car parking and picnic sites). The distance-based covariates (distances to all and major human settlements, touristic trails and forest edge) were calculated as raster proximity maps using the Euclidean distance to the nearest focal feature. The percentage of landscape (canopy) openness and average forest height were calculated based on the Digital Surface Model derived from LIDAR data (six points per m2) acquired in July 2015 within the framework of the LIFE+ ForBioSensing project (http://www.forbiosensing.pl) and provided us with a raster layer with resolution 5 m. We assumed all pixels with a vegetation height lower than 2 m were open areas and calculated their percentage share within 100 m (canopy openness at a camera trap site) and 500 m (landscape openness) grid cells. This data source was also used to describe variation in landscape topography (elevation). The percentage share of coniferous/deciduous tree species was derived through a supervised classification of a Rapid Eye satellite image acquired in June 2013 and provided by the European Space Agency (ESA; https://www.esa.int) as a raster layer with resolution 5 m. We performed the supervised classification of the Rapid Eye image using the Support Vector Machine (SVM) classifier (Mountrakis et al., 2011) as implemented in the scikit-learn Python package (Pedregosa et al., 2011). Two land cover classes were defined: coniferous dominated forest and deciduous dominated forest. All open (i.e. non-forest) areas were excluded from the classification based on the mask developed using the LIDAR data (see above). To collect spectrally homogeneous reference data for each land cover class we first ran a Large Scale Mean Shift (LSMS) segmentation as implemented in Orfeo Toolbox version 5.8.0 (Grizonnet et al., 2017). Next, the training and testing polygons were selected by image interpretation methods using very high spatial resolution data, such as the most recent available images in Google Earth (2014/10/02) and one Spot-6 scene (2015/06/27). Additionally, we used auxiliary information on tree species composition of each tree stand extracted from the State Forests inventory database (available at https://www.bdl.lasy.gov.pl). Next, we trained the SVM classifier using two-thirds of the reference dataset and the default settings in the scikit-learn implementation (Pedregosa et al., 2011). After assuring good performance of the classifier (>95% of test pixels correctly classified) we ran it for the whole study area. In the last step, we re-sampled the output raster to our focal resolutions (100 m and 500 m), calculating the percentage of coniferous and deciduous forest covering each landscape pixel. Appendix 1—figure 2 presents the correlation heatmap of all candidate spatial covariates for the ungulate models. When selecting covariates for the final models we ensured that the Pearson correlation coefficient for all pairs of included covariates was lower than 0.7. As there was a strong negative correlation (−0.97) between the percentages of coniferous and deciduous forest covering each landscape pixel we decided to include the former. For the same reason, we included distance to all settlements and excluded distance to forest edge (Appendix 1—figure 2). Also there was a strong correlation between distance to major settlements and distance to main roads (0.87), so we excluded the latter. The landscape use of lynx was determined by distance to major settlements (as indicated by the fitted model; see Figure 1 in the main text), and hence the correlation between both variables was also high (0.82). Again for the reason discussed above, we excluded the distance to major settlements from the final models.

Temporal covariates

The temporal covariates, that is temperature and snow cover, were calculated as the average values for the whole period of camera trapping at a given location. The source of the data were measurements acquired at the meteorological station in Białowieża.

Appendix 1—table 1
Summary of the two camera trapping surveys conducted during the study period.

Days – the total number of trapping days (effort), Counts – the total number of individuals recorded during independent visits (events), Trate – trapping rate. Notice the dramatic differences in the trapping rates of the carnivores between the two types of surveys. Counts of red deer female (973) and red deer male (547) do not sum to the total count (2025) as there were 505 records of individuals of undefined gender.

https://doi.org/10.7554/eLife.44937.021
SpeciesDaysCountsTrate (daily)Trate sdTrate max
Ungulates survey (05.2012–05.2014)
Wild Boar981348180.495070.8788110.16667
Red Deer981320250.202550.363403.90909
Red Deer Female98139730.096480.203912.00000
Red Deer Male98135470.055440.128671.25000
European Bison98134400.045220.232124.50000
Roe Deer98131940.020430.066110.54545
Eurasian Elk9813820.008660.047140.66667
Wolf9813510.005410.029600.36364
Eurasian Lynx981330.000270.004580.08333
Carnivores survey (09–10.2015)
Wolf30934710.191290.495346.00000
Eurasian Lynx3093350.012680.049970.33333
Appendix 1—table 2
Matrix of the foraging-related traits used for the estimation of the functional diversity index (FDis) of the entire ungulate community: body mass (based on Jedrzejewska et al., 1994) and Borowik et al., 2016), diet type (based on Hofmann, 1989); B: browsers-concentrate selectors, BG: browsers/grazers-intermediate types, G: grazers-grass roughage eaters, O: omnivores) and gut type (R: ruminants, NR: non-ruminants).
https://doi.org/10.7554/eLife.44937.022
Roe deerWild boarRed deer femaleRed deer maleMooseEuropean bison
Body mass208090150200400
Diet typeBOBGBGBG
Gut typeRNRRRRR
Appendix 1—figure 1
Heatmap of pairwise Pearson correlation values of all candidate landscape-scale covariates for the ungulates models; wolf: predicted use of the landscape by the wolf, lynx: predicted use of the landscape by the lynx, open: landscape openness, conif: % share of coniferous tree species, decid: % share of deciduous tree species, fheight: forest height, reserv: density of protected areas, elev: elevation, dsetta: distance to all settlements, dsettm: distance to major settlements, droads: distance to main roads, dedge: distance to forest edge, pois: density of tourist infrastructure.

When selecting covariates for the final models we ensured that the Pearson correlation value for all pairs of included covariates was lower than 0.7.

https://doi.org/10.7554/eLife.44937.023
Appendix 1—figure 2
The results of hierarchical clustering on the principal component analysis.

The five clusters were identified by assessing the inertia (i.e. change in within cluster homogeneity) gained by cutting the tree at different levels and the ecological interpretability of the resulting clusters. The colors of the clusters correspond to those in Figure 10 in the main text.

https://doi.org/10.7554/eLife.44937.024
Appendix 1—figure 3
Spatial patterns of human activity in Białowieża Forest derived from the publicly available world-wide data provided by STRAVA (https://www.strava.com/heatmap) combined with (A) predicted use of the landscape by the wolf, (B) predicted use of the landscape by the lynx, (C) predicted use of landscape by the red deer (both sexes) and (D) distance to major settlements used in the models for both large carnivores.

The dark-red thick lines indicate high human activity and the light-blue thin lines indicate low human activity. For (A), (B) and (C) the color gradient from blue to red corresponds to the gradient from low to high use of the landscape by the presented species.

https://doi.org/10.7554/eLife.44937.025
Appendix 1—figure 4
Comparison of predictions of wolf space use with existing radio-tracking data from collared wolves in Białowieża Forest (period 1994–1999, Jędrzejewski et al., 2007) and observations of known den locations (period 1993–2007, MRI PAS unpublished data).

In the background we plotted the raster map of wolf space use predicted by our model together with the raw camera trapping data (‘bubbles’ with trapping rates at each location). On top of it we plotted two datasets: 1) the locations of the core areas (dashed red circles) of the annual territories of the four wolf packs present in the area (50% MCP based on telemetry) and 2) den locations.

https://doi.org/10.7554/eLife.44937.026
Appendix 1—figure 5
Bubble plot presenting raw camera trapping data for lynx collected during the carnivore survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.027
Appendix 1—figure 6
Bubble plot presenting raw camera trapping data for red deer female collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.028
Appendix 1—figure 7
Bubble plot presenting raw camera trapping data for red deer male collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.029
Appendix 1—figure 8
Bubble plot presenting raw camera trapping data for wild boar collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.030
Appendix 1—figure 9
Bubble plot presenting raw camera trapping data for roe deer collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.031
Appendix 1—figure 10
Bubble plot presenting raw camera trapping data for bison collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.032
Appendix 1—figure 11
Bubble plot presenting raw camera trapping data for moose collected during the ungulate survey.

Each bubble is a camera trap location and its size is proportional to the daily trapping rate (i.e. no. of individuals observed/no. of days).

https://doi.org/10.7554/eLife.44937.033
Appendix 1—figure 12
Basic meteorological data for the study period (05/2012 – 09/2014): 1) averaged daily temperature [C], 2) precipitation [mm] and 3) snow cover [cm].
https://doi.org/10.7554/eLife.44937.034
Appendix 1—figure 13
Distribution of camera trap locations grouped by quality of view in front of a camera and different UNESCO zones in Białowieża Forest (1: Strict Reserve of BNP, 2: BNP and nature reserves, 3: valuable unprotected forest stands, 4: managed forest stands).

Each camera location was manually tagged with a 4-level categorical label (‘Exclude’, ‘Acceptable’, ‘Good’, ‘Perfect’) describing the quality of view in front of the camera. Locations labeled with the ‘Exclude’ were not included in further analysis and are not shown on this figure.

https://doi.org/10.7554/eLife.44937.035
Appendix 1—figure 14
Red deer female – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.036
Appendix 1—figure 15
Red deer female – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.037
Appendix 1—figure 16
Red deer male – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.038
Appendix 1—figure 17
Red deer male – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.039
Appendix 1—figure 18
Roe deer – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.040
Appendix 1—figure 19
Roe deer – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.041
Appendix 1—figure 20
Wild boar – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.042
Appendix 1—figure 21
Wild boar – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.043
Appendix 1—figure 22
European bison – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.044
Appendix 1—figure 23
European bison – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.045
Appendix 1—figure 24
Moose – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.046
Appendix 1—figure 25
Moose – detection rate; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.047
Appendix 1—figure 26
Wolf – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.048
Appendix 1—figure 27
Eurasian lynx – landscape use; the posterior distributions of all parameters.
https://doi.org/10.7554/eLife.44937.049
Appendix 1—figure 28
Red deer female – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.050
Appendix 1—figure 29
Red deer male – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.051
Appendix 1—figure 30
Roe deer – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.052
Appendix 1—figure 31
Wild boar – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.053
Appendix 1—figure 32
European bison – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.054
Appendix 1—figure 33
Moose – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.055
Appendix 1—figure 34
Wolf – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.056
Appendix 1—figure 35
Eurasian lynx – model diagnostic plots.
https://doi.org/10.7554/eLife.44937.057
Appendix 1—figure 36
Estimated effects of the covariates from the model for all ungulate species using only a subset of camera trap data covering a three month period (August - October) closely matching the period of the carnivore survey.

The covariates explain the spatial variation in the parameter λ, which is the expected number of individuals using a given landscape grid cell (25 ha pixel) during the sampling period (i.e. the relative density; see the general and formal description of the model in the Materials and methods section).

https://doi.org/10.7554/eLife.44937.058
Appendix 1—figure 37
The spatial predictions for the parameter λ, which is the expected number of individuals using a given landscape grid cell (25 ha pixel) during the sampling period (see the general and formal description of the model in the Materials and methods section).

The predictions are based on the model for all ungulate species using only a subset of the camera trap data covering a three month period (August - October) closely matching the period of the carnivore survey.

https://doi.org/10.7554/eLife.44937.059
https://doi.org/10.7554/eLife.44937.020

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
    R: A Language and Environment for Statistical Computing
    1. R Core Team
    (2019)
    Vienna, Austria: R Foundation for Statistical Computing.
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
    Vegetation Dynamics in Temperate Lowland Primeval Forests: Ecological Studies in Białowieża Forest, Geobotany Series
    1. J Faliński
    2. K Falińska
    (1986)
     Dordrecht: Springer.
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
    Data Analysis Using Regression and multilevel/Hierarchical Models
    1. A Gelman
    2. J Hill
    (2006)
    Cambridge University Press.
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
    Predation in Vertebrate Communities, Ecological Studies
    1. B Jędrzejewska
    2. W Jędrzejewski
    (1998)
    Berlin, Heidelberg: Springer.
  49. 49
    Large carnivores and ungulates in European temperate forest ecosystems: bottom-up and top-down control
    1. B Jędrzejewska
    2. W Jędrzejewski
    (2005)
    Large Carnivores and the Conservation of Biodiversity pp. 230–246.
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
    Applied Hierarchical Modeling in Ecology: Analysis of Distribution, Abundance and Species Richness in R and BUGS: Volume 1:Prelude and Static Models
    1. M Kery
    2. JA Royle
    (2016)
    Academic Press.
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
    Vegetation landscapes of bialowieza forest
    1. W Kwiatkowski
    (1994)
    Phytocoenosis 6:35–88.
  72. 72
  73. 73
  74. 74
  75. 75
    FactoMineR : an R package for multivariate analysis
    1. S Lê
    2. J Josse
    3. F Husson
    (2008)
    Journal of Statistical Software, 25, 10.18637/jss.v025.i01.
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
    Scikit-learn: machine learning in Python
    1. F Pedregosa
    2. G Varoquaux
    3. A Gramfort
    4. V Michel
    5. B Thirion
    6. O Grisel
    7. M Blondel
    8. P Prettenhofer
    9. R Weiss
    (2011)
    Journal of Machine Learning Research 12:2825–2830.
  93. 93
  94. 94
  95. 95
  96. 96
    QGIS geographic information system
    1. QGIS Development Team
    (2017)
    Open Source Geospatial Foundation.
  97. 97
  98. 98
  99. 99
  100. 100
  101. 101
  102. 102
  103. 103
  104. 104
  105. 105
  106. 106
  107. 107
  108. 108
  109. 109
  110. 110
  111. 111
  112. 112
  113. 113
  114. 114
  115. 115
  116. 116
  117. 117
  118. 118
  119. 119
  120. 120
  121. 121
  122. 122
    The niche, biogeography and species interactions
    1. JJ Wiens
    (2011)
    Philosophical Transactions of the Royal Society B: Biological Sciences 366:2336–2350.
    https://doi.org/10.1098/rstb.2011.0059
  123. 123
  124. 124
    Global climate and the distribution of plant biomes
    1. FI Woodward
    2. MR Lomas
    3. CK Kelly
    (2004)
    Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 359:1465–1476.
    https://doi.org/10.1098/rstb.2004.1525
  125. 125
  126. 126
  127. 127

Decision letter

  1. Ian T Baldwin
    Senior Editor; Max Planck Institute for Chemical Ecology, Germany
  2. Christian Rutz
    Reviewing Editor; University of St Andrews, United Kingdom
  3. Christian Rutz
    Reviewer; University of St Andrews, United Kingdom
  4. Joern Theuerkauf
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Thank you for submitting your article "Linking spatial patterns to trophic interactions in terrestrial herbivore communities" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Christian Rutz as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Ian Baldwin as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Joern Theuerkauf (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This study potentially makes a valuable contribution to our understanding of the processes structuring herbivore communities. While generally a well-studied topic, the study provides rare data from a temperate forest system, identifies functionally-distinct 'herbiscapes' across the landscape, and attempts to measure the influence of human activity. That said, the reviewers have raised a relatively large number of major concerns that need to be carefully addressed in a revision.

Essential revisions:

1) Data coverage. One of the main concerns is that the datasets produced by your camera-trapping approach (as summarised in Appendix 1—table 1) may be too small to enable robust statistical inferences. Counts are relatively small for all but two of the ungulate species (wild boar and red deer), and only few detections were made for the two carnivores (even during the repeat survey). In fact, the data for lynx were so sparse, with just 35 detections, that we had little confidence in the relationships shown in Figure 3D-F, and we actually remained unconvinced that this species should be considered a major predator of the ungulates studied (except, perhaps, for their young, but this would require further justification and analyses). While a better sample was achieved for wolves, with 471 detections, we were concerned that this may still be insufficient to map the species' ranging behaviour. In fact, we were surprised that you had not used data from extensive earlier radio-tracking projects in the study area. To address these concerns, please: (a) remove the lynx from your species set (toning down claims accordingly, about this being a 'multi-predator' system); (b) provide further analyses demonstrating that detection rates for all remaining species were sufficient to allow robust statistical inferences; and (c) explicitly compare your camera-trap data for wolves to earlier radio-tracking data, and assess whether using the latter yields similar conclusions overall.

2) Earlier work. We were surprised to find no mention of a previous study in the Bialowieza Forest, which was very similar in scope (Theuerkauf J and Rouys S 2008, Habitat selection by ungulates in relation to predation risk by wolves and humans in the Bialowieza Forest, Poland. Forest Ecology and Management 256, 1325-1332). Please provide a detailed comparison of your results to those of this earlier study (which had combined pellet counts on transects, with radio-tracking data). It would be reassuring if both studies independently reached similar conclusions, based on different methodology, and if there are disagreements, these should be highlighted and explored further.

3) Camera-trapping methodology. Detection rates are expected to be highly dependent on the placement of camera traps, which can affect inferences about species' habitat use. The area covered by the camera is usually larger in habitats with low ground and scrub cover (not the canopy cover; remote sensing can therefore not detect these differences), which can bias habitat-selection estimates towards areas with high visibility. For example, looking at the distribution map for red deer (Appendix—figure 3), we were surprised that the highest densities were not observed in the strict reserve. Since ground vegetation cover seems higher in the strict reserve than the commercial forest, as noted by one of the reviewers, the camera-trapping method could have biased the distribution of red deer towards the commercial forest. There have been extensive radio-tracking studies on wolves, lynx, bison and red deer in the Bialowieza Forest, and these should be consulted for formal, quantitative comparison.

4) Correlative evidence. It is important to acknowledge that your study draws inferences from correlative evidence, while establishing causality would require experimental manipulation of key parameters. Throughout the manuscript, terms are used that imply causality (to drive, determine, affect etc.), and these should be replaced with appropriate alternatives (to correlate, be associated with etc.) or combined with suitable qualifiers (may, seems, appears to be etc.). Two analyses in particular require more nuanced wording. First, the 'herbiscape' analysis is identified as a primary point of novelty, and while conceptually this is correct, the evidence provided is largely correlative in nature. In the absence of direct experimental evidence of the effects of herbivores on the specific plant species examined (through controlled exclosure experiments), it is impossible to conclude that the spatial distribution of herbivores is causing the observed vegetation patterns. We recognise that such evidence is difficult to obtain, and we are not suggesting that additional work is required, but we think that this issue needs to be acknowledged (e.g., in subsection “Spatial heterogeneity in the landscape distribution of large herbivores” of the Discussion). Second, while it is clear from the data that red deer did not use the same areas as wolves, it is not possible to infer that they 'avoided' those areas (which implies causality). In fact, wolves may only select areas of high red deer density during periods of hunting; the negative correlation between deer and wolf distributions could therefore arise if areas used by wolves when they are not hunting are not preferred habitats of deer. Besides, as red deer are the major prey of wolves in the study area, why would wolves then not select areas of higher red deer density? We appreciate that experimental manipulation of this system is impossible, but note that analyses of spatio-temporal patterns (e.g., from simultaneous radio- or GPS-tracking of wolves and deer) would have considerably strengthened inferences. Please tone down your discussion of wolf-deer relationships accordingly.

5) Presentation. This is a complex study, both in terms of the biological processes studied and the analysis techniques used, and it is critical to help readers follow the narrative. This is especially important as eLife's readership is very broad, and papers should be intelligible to non-specialists. To improve the presentation of the study, we suggest three main revisions. (a) Technical terms should be avoided where possible, or at least clearly defined at the outset. For example, the Introduction talks about "cascading effects", "hyper-keystone species" etc. While these may be frequently-used terms in this research field, without further explanation, many non-specialist readers may struggle to understand the intended meaning. In general, the Introduction could be more focussed and reader-friendly. (b) Since the Results and Discussion sections precede the Materials and methods, some aspects of the study system and methodology remain unclear until much later in the narrative. For example, the species under investigation haven't been mentioned until they are referred to in various parts of the Results section, and 'herbiscapes' are discussed in the Results before the concept is explained in the Discussion. There are two possibilities for dealing with this. Either the Materials and methods are moved forward, for a traditional I-M-R-D structure, or the front end of the narrative is carefully revised, to equip the reader with all the information they need to understand the results and their interpretation. (c) We think it would really help readers if there was a 'graphical abstract' that illustrates, in two separate panels: (i) how data were collected and processed; and (ii) the biological relationships uncovered by the analyses (humans --> wolves --> prey --> vegetation). Panel (ii) could provide explicit call-outs to figures and materials in the online supplement that provide support for a given link, helping readers navigate these various materials.

6) Human activity. Did your camera traps record human activity (on/off established tracks)? If so, why were these data not included in the models, instead of (or in addition to) your chosen proxy variable (distance to settlements)? Appendix—figure 3 shows human activity data, but it is unclear how these were obtained, and "combined" with other datasets, and whether any of this was used in the models. Incidentally, human detections should enable a very useful assessment of the reliability of your camera-trapping method (the distribution of human observations is expected to match the known distribution of human activity in the forest); see also comment #3 above.

7) Sampling periods. The analyses assessing patterns for a subsampled ungulate dataset (for the months covered by the carnivore survey) are very useful. But, it is not clear why an extra month was included here (Aug-Oct; subsection “Carnivores survey”). Why not use Sept-Oct, to make it strictly comparable? More importantly, you say that the "results were similar to those based on the full dataset", yet careful comparison between Figures 4/5 and Appendix 1—figures 36/37 reveals important differences. This should be acknowledged, further investigated, and discussed. Which of your main conclusions remain unchanged when you use this subsampled dataset, and which ones are affected?

8) Broader implications. The study makes a notable contribution to our understanding of key processes in community ecology, but we felt it would be useful to discuss some broader implications in the Conclusions section. Specifically, are any of the findings relevant for ongoing conservation/management efforts, in the study area or elsewhere?

9) Transparent reporting. Some of the points raised in eLife's guidance document require attention (e.g., power analyses, which relates to comment #1 above).

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Linking spatial patterns of terrestrial herbivore community structure to trophic interactions" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by Christian Rutz as the Reviewing Editor and Ian Baldwin as the Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Joern Theuerkauf (Reviewer #3).

The reviewers have carefully evaluated your revised manuscript and response letter. While reviewer #2 was happy with the revisions overall, apart from a few minor points that still need addressing, reviewer #3 remained concerned about data quality, and specifically, the mismatch in timing of the predator and prey datasets. The reviewers' original reports are appended below for reference. We share reviewer #3's concerns, but are prepared to consider a final revision that addresses these points more explicitly and tones down claims further. It would also be important to better acknowledge throughout that correlations between spatial species distributions were examined, not interactions (as for example implied by the use of the term “cascading effects”).

Reviewer #2:

All of my major and minor comments have been addressed in this revised version of the manuscript. I believe the revision is much-improved and would recommend publication after addressing the few minor comments below.

Minor Comments:

Introduction paragraph two: Unsure of the accuracy of this claim that this subject has often been ignored – "The strong spatial variation in herbivory impact, resulting from the space use of by different functional groups of herbivores has often been ignored." Suggest re-wording to simply state the effect, as it has been shown elsewhere, e.g. "There is often strong spatial variation in herbivory impact resulting from the space use of by different functional groups of herbivores."

Subsection “Spatial heterogeneity in the landscape distribution of large herbivores” paragraph three: Avoid repeated and strong claims of novelty (as this type of study has been done in other systems, as you say). Suggest re-wording to "Similar studies in African ecosystems have shown that…"

In the same subsection: Suggest rewording to focus specifically on the benefits of this study from providing information on a temperate ecosystem – this is the major point of novelty and, in my view, what makes this study a valuable contribution to this body of literature, and should be emphasized. The current wording does not make this clear.

Reviewer #3:

The first concern regarding the manuscript was the robustness of statistical inference. In the rebuttal, the authors assume that if there is no problem with the robustness of the statistical analyses, then the results are reliable. However, the problem does not lie with the analyses but with the quality of data. If the data are not representative in the first place, then there is even no way to verify the robustness of the results. The major problem I see is the temporal mismatch of data sampling. While you sampled ungulate distribution over 2 years, you only sampled predators during a single period of 2 months. Besides, you state that "In August-September, pups begin to travel with other pack members and move more widely through their territory, returning to the core of their territory on a regular basis (Jędrzejewski et al., 2001)". This is however not completely correct. Wolves actually use core areas only during the pup rearing season but not during the rest of the year. Therefore, seasonal data is not representative for the year-round space use. You can therefore compare space use by wolves in September-October only with prey space use during the same season. Otherwise, it is possible you will get random results. Even if prey space use in September-October is representative for the whole year (using the whole year will actually create only the impression of more robust data by artificially increasing sample size, but as the predator sampling is reduced to a short period, it does not help), wolf space use in September-October is definitely not representative for the whole year and therefore also not the best period to quantify space use as you state at the end of the first paragraph in subsection “Carnivores survey”. I would therefore not agree with your statement in the rebuttal regarding seasonal differences in wolf space use: "We appreciate this idea but we did not have data on large carnivores from other seasons. Moreover, we expect little variation in wolf use of the landscape". On the contrary, there are huge differences in wolf space use comparing spring/summer and autumn/winter. This is especially important as September-October also overlap with the rut of red deer, thus the period when wolves rather select hunting male and not female deer (during this period you found that females more correlated with wolves). Given that you did also not consider any temporal information (times of hunting vs times of resting), it is not really possible to assess if the correlations are actually based on real predator-prey interactions or if they are random results. It is possible that the results actually are results of trophic relations, but the data structure does not allow explicit conclusions.

https://doi.org/10.7554/eLife.44937.061

Author response

Essential revisions:

1) Data coverage. One of the main concerns is that the datasets produced by your camera-trapping approach (as summarised in Appendix 1—table 1) may be too small to enable robust statistical inferences. Counts are relatively small for all but two of the ungulate species (wild boar and red deer), and only few detections were made for the two carnivores (even during the repeat survey). In fact, the data for lynx were so sparse, with just 35 detections, that we had little confidence in the relationships shown in Figure 3D-F, and we actually remained unconvinced that this species should be considered a major predator of the ungulates studied (except, perhaps, for their young, but this would require further justification and analyses). While a better sample was achieved for wolves, with 471 detections, we were concerned that this may still be insufficient to map the species' ranging behaviour. In fact, we were surprised that you had not used data from extensive earlier radio-tracking projects in the study area.

We appreciate the concerns of the reviewers regarding the robustness of the statistical inferences in relation to the sample size. However, this comment is too broad to address here, as the reviewer had not identified any specific technical problems with our statistical approach (relating to, e.g., a particular part of the methods, the statistical framework used, formal conception of the model).We are convinced that we can defend the structure and specification of our model and its fit to the data. We provided a very detailed description of the applied statistical framework, including a general and formal description of the model, its specification and implementation details together with the data and the source code of our analysis in Python. Our model was built on a well-established statistical framework (Guillera-Arroita et al., 2011, Guillera-Arroita et al., 2012; Hughes and Haran, 2012; Johnson et al., 2013; Kery and Royle, 2016; MacKenzie et al., 2002; Royle et al., 2007; Royle and Nichols, 2003). In the supplementary materials we have provided the standard diagnostic plots for all models of all species (both carnivores and ungulates) which show, without doubt, that our models, including the lynx model, do fit to the data (see Appendix 1—figures 28-35). In the supplementary materials an interested reader can also find the figures presenting the posterior distributions for all estimated parameters (for each species; Appendix 1—figures 14-27) as well as track-plots with the MCMC chains used for sampling from the posterior distributions. Therefore, we are not sure about the validity of the reviewer’s criticism on this point.

Species with a low count, such as the lynx, logically have a larger uncertainty in their parameter estimates as expressed by the larger size of the credible intervals and overlap with zero. Importantly, we want to strongly emphasize that in addition to the camera trap records of animals, all the recorded zeros are equally important. Our model takes all this information (both zero and non-zero counts) into account. A zero output can result from an ecological process or can be caused by imperfect detection. Our hierarchical model takes both these possibilities into account (see Materials and methods subsection “Formal description”). One model layer explains the variation in the data arising from ecological processes (spatial variation in density of animals), and another layer deals with both ecological (e.g. habitat selection, movement, seasonal activity levels) and observational processes (e.g. detection issues, camera failures). Additionally, we find it worth emphasizing once more that we used data from 894 locations resulting in 9813 trap nights from the ungulate survey, and from 73 locations resulting in 3093 trap nights from the carnivore survey (see Appendix 1—table 1). This should have resulted in robust estimates of animal space use and detection rates; see for example Parsons et al., 2018.

Taking all the above into account, we are convinced that our data and statistical approach allows for a robust statistical inference, and we therefore see no reason how it could be “insufficient to map the species’ ranging behaviour”. Below we present more specific responses to the three points raised by the reviewers.

To address these concerns, please: a) remove the lynx from your species set (toning down claims accordingly, about this being a 'multi-predator' system);

We agree that the sample size for the lynx is particularly small, but we do not see a reason to exclude it from our analysis, firstly because of the technical issues discussed above, and secondly, and even more importantly, for ecological reasons. That the Eurasian lynx is a relevant predator in the studied system has been shown by several studies (Jędrzejewska and Jędrzejewski, 1998; Okarma et al., 1997). These studies showed that the lynx’s major prey is the roe deer, but that the red deer also constitutes between 22 and 61% of its diet (Okarma et al. 1997, Jędrzejewski et al., 1993). The annual mortality caused by lynx constitutes 21-36% of the spring numbers of roe deer and 6-13% of the spring number of red deer (Okarma et al., 1997). Lynx predation is therefore a major mortality factor for roe deer and a secondary factor for red deer (Okarma et al., 1997, Jędrzejewski et al., 1993, Jedrzejewska et al., 1997, Jędrzejewska and Jędrzejewski, 2005). Therefore, the evidence fully supports the lynx being an important predator in this system (as we discussed in paragraph six of subsection “The effect of large carnivores on the spatial structure of the large herbivore community”), despite their low numbers. Moreover, our experimental studies have shown that red deer clearly react with anti-predatory behavior to olfactory cues of the lynx in the Białowieża Forest (Wikenros et al., 2015). As the present study shows, these effects on their prey species do not create observable spatial patterns at the landscape scale (as in the case of the wolf). This lack of a clear effect on the landscape use of its prey species likely results from the combined effects of lynx’s low density and low site fidelity (Schmidt et al., 1997. Poland. Acta Theriologica 42:289–312. DOI: 10.4098/AT.arch.97-30) and the previous observation that lynx activity is explained by even finer-scale habitat characteristics than those studied in our camera trapping survey (Podgórski et al., 2008).

We added a new paragraph in the Discussion emphasizing that lynx is a relevant predator in our system and discussing the lack of clear effect on the spatial pattern of landscape use by the prey species of lynx (paragraph six of subsection “The effect of large carnivores on the spatial structure of the large herbivore community”).

b) provide further analyses demonstrating that detection rates for all remaining species were sufficient to allow robust statistical inferences; and

We have no doubts that our statistical inference was robust. Please, see our general response above where we elaborate on this issue. However, if the reviewers are not fully convinced by our arguments we would like to refer them to a recently published paper by Parsons et al., 2018, in eLife: “Our sample size goal was 20 spatial replicates (equating to ~420 trap nights), which has been found to maximize precision for estimating detection rate (Kays et al., 2010; Rowcliffe et al., 2008)”. In comparison, in our study we had 894 spatial replicates resulting in 9813 trap nights for the ungulate survey (see Appendix 1—table 1). Therefore, our study is far more data-rich and we feel that the reviewers’ lack of confidence in our study is out of place in this context, especially because no formal part of our analysis was explicitly questioned.

c) explicitly compare your camera-trap data for wolves to earlier radio-tracking data, and assess whether using the latter yields similar conclusions overall.

The existing radio-tracking data from collared wolves in Białowieża Forest was collected over 20 years ago (1994-1999) and since then many things have changed in the forest. For example, the enlargement of Białowieża National Park in 1996, changes in forest management practices, the development of tourism and its infrastructure in the region and last, but not least, a mass dieback of spruce stands in recent years (Boczoń et al., 2018; Mikusiński et al., 2018). As wolves base their habitat selection on human activity, habitat characteristics and prey availability, the changes that have occurred since the end of the 1990s have likely affected wolf space use. Moreover, core areas of wolf territories (i.e. where they reproduce) are not fixed in space and time because female wolves can use different dens each year, and den locations can even change by several kilometers within a season (Schmidt et al., 2008). Therefore, the position of the current wolf core areas could have significantly changed from the telemetry studies observations from the late 1990s (see also Kuijper et al., 2015).However, previous work in Białowieża Forest has showed that human related factors are a key factor determining wolf space use (Theuerkauf et al., 2003) and we know that these factors have been present in a similar spatial arrangement for decades or even centuries (see subsection “Specification - wolf and lynx model” in the Materials and methods section and paragraph seven in subsection “The differential effect of large carnivores on the spatial structure of the large herbivore community” of our Discussion); thus, we expect a general agreement between our camera trapping and the radio-tracking data. To check this, we compiled a new map (Appendix 1—figure 14). In the background we plotted the raster map of wolf space use predicted by our model together with the raw camera trapping data. On top of it we plotted two datasets: 1) the locations of the core areas of the annual territories of four wolf packs based on telemetry studies in the 1990s (Jędrzejewski et al., 2007) and 2) observations of known den locations (period 1993–2007) and wolf howling during the reproductive season (2000–2014) (Nowak et al, 2007 and unpublished data). The map shows that there is a good spatial agreement between the previous telemetry data and the spatial predictions of our model. The only difference is the activity center of the NW wolf pack, which according to the camera trapping data has moved to the north part of Białowieża Forest when compared to the radio-tracking data. This is in line with the most recent (available) locations of dens in this area, which also occurred further north. In summary, our model fitted to the camera trapping dataset (471 wolf detections) is highly reliable as it not only conforms to the general pattern of the wolves’ space use determined by telemetry 20 years ago, but also likely reflects the wolves’ response to ongoing environmental changes in the study area. In other words, using the outdated telemetry data for this study would not be justified.

We compiled a new figure (Appendix 1-figure 4) and added a short paragraph referring to this good agreement in the main text (paragraph two in subsection “Large carnivores”). Moreover, we placed the entire text of this particular response in the appendix as we believe it is a good description of the rationale behind this comparison of our results and the radio-tracking data from ‘90s.

2) Earlier work. We were surprised to find no mention of a previous study in the Bialowieza Forest, which was very similar in scope (Theuerkauf J and Rouys S 2008, Habitat selection by ungulates in relation to predation risk by wolves and humans in the Bialowieza Forest, Poland. Forest Ecology and Management 256, 1325-1332). Please provide a detailed comparison of your results to those of this earlier study (which had combined pellet counts on transects, with radio-tracking data). It would be reassuring if both studies independently reached similar conclusions, based on different methodology, and if there are disagreements, these should be highlighted and explored further.

We have to admit that we have indeed not carefully enough considered this work, mainly due to its completely different methodology. However, this reviewer’s comment has stimulated us to re-examine the paper by Theuerkauf and Rouys, 2008, and compare our research approaches to further validate our results. We have now included a short paragraph in which we refer to this paper and discuss the additional insights that our study has yielded, see lines paragraph five of subsection “The differential effect of large carnivores on the spatial structure of the large herbivore community”. Theuerkauf and Rouys, 2008, generally found that:

1) Habitat alteration by forest exploitation and human hunting influenced the density distribution of ungulates more than predation risk by wolves.

We basically showed similar patterns, i.e. that humans shape the ecological interactions in Białowieża Forest and that these cascade down trophic levels. However, we explored these patterns in more detail, decomposing the inter-specific variation in ungulate space use by using fine-scale, spatially explicit explanatory layers. We were able to detect species specific relationships that included human-predator-ungulate interactions. Moreover, we extrapolated this to the theoretically robust ecological model of herbiscapes.

2) They also showed that (wolf) prey density was not higher in different wolf risk zones, regardless of human hunting.

We showed that the red-deer’s (the main prey of wolf) use of landscape was spatially negatively associated wolf’s high-use areas. As Theuerkauf and Rouys concluded, the “predation risk becomes more complex in multi-predator systems”, our study provided much deeper insight into this complexity showing an interactive, context-dependent character of these relationships.

In general, our results correspond to the findings of Theuerkauf and Rouys, 2008, at least in the comparable parts. However, our study was based on a much finer scale inference and different methodology and therefore there are also differences between the studies. We focused our inference on the variation in ungulate distributions over the landscape, which is a spatial process manifesting itself in observable spatial patterns emerging from interacting ecological mechanisms operating at a range of spatial scales (Levin, 1992). To gain insight into this process we recorded and reconstructed these spatial patterns by sampling randomly with respect to species and the hypothetical mechanisms driving their distributions (habitat, humans, wolf). We then tested the importance of these mechanisms via a regression-like analysis. We used an intensive large-scale camera trap study to collect detailed, spatially-explicit information on species distribution. Further, we combined these observations with GIS and high-resolution remote sensing data and used a spatially-explicit hierarchical modeling approach to explain the variation in the collected camera trapping data. As a result, our study greatly expand upon the findings of Theuerkauf and Rouys, 2008, study by adding extra resolution both in space and in species level effects. We gained additional insights due to the higher resolution of the camera trap data we used and a (close to) continuous landscape classification based on satellite images and available GIS data.

Moreover, some of the differences between Theuerkauf and Rouys and our results could be attributed to the fact that we were able to account for the contemporary, highly heterogeneous conservation zoning in BF, which has increased significantly since their work. As a result, the areas that were previously classified as “exploited” forest, in reality constitute well preserved habitats that are now excluded from exploitation.

3) Camera-trapping methodology. Detection rates are expected to be highly dependent on the placement of camera traps, which can affect inferences about species' habitat use. The area covered by the camera is usually larger in habitats with low ground and scrub cover (not the canopy cover; remote sensing can therefore not detect these differences), which can bias habitat-selection estimates towards areas with high visibility.

We agree with the reviewer that detection rates depend on many factors, including near-ground visibility in front of a camera. This could potentially have been an issue in the ungulate survey, where camera traps were placed pseudo-randomly within the forest. For this reason, to standardize recording conditions at each camera trap site, we always tried to mount a camera assuring a clear view of at least 20 m, as written in the Materials and methods section. Whenever this was impossible, we randomly looked for an acceptable place to mount a camera as close as possible to the location given by the coordinates pre-computed prior to the field work in QGIS. We have now added this clarification to the Materials and methods section. In our model we specifically took into account factors that could have affected detection probability. As such, our detection rate sub-model was designed to control for other sources of variation in the detectability of animals relating to both ecological (e.g. habitat selection, movement, seasonal activity levels) and observational processes. Please, see the Materials and methods section for more details. Additionally, in our camera trapping dataset, each camera location was tagged with a 4-level categorical label (“Exclude”, “Acceptable”, “Good”, “Perfect”) describing the quality of view in front of the camera. We used this data to prepare a new figure showing that there is no difference between the managed part of BF and BNP (Appendix 1—figure 13).

For example, looking at the distribution map for red deer (Appendix 1—figure 3), we were surprised that the highest densities were not observed in the strict reserve. Since ground vegetation cover seems higher in the strict reserve than the commercial forest, as noted by one of the reviewers, the camera-trapping method could have biased the distribution of red deer towards the commercial forest.

We would like to emphasize once more that the strict reserve and the commercial forest are not homogeneous and discrete forest patches. There is large variation in habitat structure in both areas as illustrated e.g. by our high resolution satellite data (see Figure 2). The commercial forest is also a patchwork of (often extremely) variable forest stands characterized by diverse abiotic conditions (Kwiatkowski, 1994) and forest management history. Despite intensive timber harvest in the past, many forest patches remain natural or of semi-natural character. The large variation in forest canopy cover (see Figure 2) also results in large variation in ground vegetation cover. Consequently, the density of the understory is equal or often greater in certain areas of the commercial forest than in the strict reserve.

Concerning deer densities; the southern part of the Strict Reserve was actually one of the many red deer density hotspots when considering the entire Białowieża Forest (Figure 5, 6, 9). This interesting pattern of higher densities of red deer in the south part of the strict reserve, outside the wolf core area, is in line with our previous studies based on red deer visitation rates and behaviour (Kuijper et al., 2015), and spatial patterns of browsing intensity (Kuijper et al., 2013), which also showed these clear gradients in deer abundance.

There have been extensive radio-tracking studies on wolves, lynx, bison and red deer in the Bialowieza Forest, and these should be consulted for formal, quantitative comparison.

We have noticed that a large part of the reviewers’ critique of our work is related to the quality of our camera trapping data (points #1, #3, #6) and lack of a formal and quantitative comparison of our results to telemetry-based studies (points #1, #3, #4), and so would like to briefly discuss these issues.

Camera trapping, although not a panacea for surveying wildlife (Wearn and Glover-Kapfer, 2019), is a well-established, non-invasive method of collecting field data on animal abundance, distribution, temporal activity and space use (Burton et al., 2015; Karanth et al., 2017; Pfeffer et al., 2017; Wearn and Glover-Kapfer, 2017). Camera trapping provides reliable data for population- and community-level inferences by itself and there is no a priori requirement to complement camera trapping with telemetry data, although the latter can provide important auxiliary information e.g. for mark-recapture studies (Dillon and Kelly, 2008; Soisalo and Cavalcanti, 2006).

We fully recognize the value of telemetry, especially GPS-based, which is a powerful method providing rich information on animal movement ecology, activity patterns, resource selection etc. However, it is also not a panacea for ecological studies as pointed out by Hebblewhite and Haydon, 2010, in their critical review. One of their conclusions was that one of major disadvantages of GPS-telemetry are the generally small sample sizes and poor population-level inferences. For studies of resource selection, sample size requirements of more than 30 units for robust population-level inferences likely still apply, but still this recommendation does not address representativeness of the sample, which is a function of the total size of the population (Hebblewhite and Haydon, 2010). Therefore, radio-tracking studies from Białowieża forest from 12 collared wolves (Jędrzejewski et al., 2007) out of a wolf population of ~30 (4 packs 6-7 members each) is not the same as data from 19 collared red deer (Kamler et al., 2007) out of a red deer population of ~3000 (5-6 individuals/km2 x 580km2) or 29 collared wild boars (Podgórski et al., 2013) out of a wild boar population of ~3000 (5-6 individuals/km2 x 580km2). Besides this generally low number of collared individuals as a proportion of the entire population (in the case of ungulates), spatial coverage is also essential. The collared individuals of both red deer and wild boar were captured and radio-tracked in a limited area of the Białowieża Forest (at the border of National Park). So the observed patterns from telemetry studies are representative for only a very small part of the population from a particular part of the Białowieża Forest landscape.

The main ecological process we studied, the distribution of five ungulate species over the BF landscape, is a spatial process. Our sampling scheme was generally designed to maximize the probability of capturing the overall landscape scale spatial variation in local abundances of each species. We recorded and reconstructed spatial patterns of all the study species by sampling the continuous BF landscape with camera traps placed randomly with respect to species and the hypothetical mechanisms driving their distributions (habitat structure, humans, wolf). That is, we used an intensive, large-scale and high-resolution camera trap network covering the entire study area to collect detailed, spatially-explicit information on species distributions. We argue that this kind of study is impossible to do with telemetry. To record the targeted spatial patterns, observations had to be made at the level of the entire population of each species and not at the level of single individuals. No matter how accurately tracked in space and time, data from GPS-collared individuals would not allow the re-constructing of approximated density-surfaces, which are the products of landscape use of the entire population of a species. The only exception from this is the radio-tracking data for wolves. In this case the recorded landscape use by collared individuals (c. 50% of all wolves in BF, with each pack represented) is a very good approximation of the landscape use by the entire wolf population. These are the most fundamental reasons why we could not use telemetry and why we cannot formally nor quantitatively compare our results with radio-tracking data from the past, except in the case of wolves (but please see our answer to #1). Finally, our camera trapping and analytical approach, as explained above, enabled us to obtain a reliable picture of wolf spatial patterns, comparable to previous telemetry results, while simultaneously providing high resolution data on all species of the large mammal community with the same unified methodology. We have clarified this in the Materials and methods section.

However, what we can do is a partial comparison of our results (i.e. from factors influencing detection rates of single species to patterns observed at the community level) with the previous studies from Białowieża forest, including those based on telemetry. The idea behind this is similar to the philosophy behind pattern-oriented modelling, where multiple patterns observed in complex systems at different hierarchical levels and scales are used systematically to optimize model complexity and reduce uncertainty (Grimm, 2005). We believe we did a good job of this in the results and Discussion sections in the main text. Here we only briefly list some examples of the partial results and patterns observed in our data that are directly comparable to the previous studies from Białowieża Forest, and we refer to the figures presenting these results:

1. The use of the landscape by wolves is primarily determined by human related factors, as was shown earlier by Theuerkauf et al. (2003); Figure 3,

2. The distribution of red deer, the main prey of the wolf (Jędrzejewski et al., 2002), negatively correlates with the landscape use of wolves, as was shown by Kuijper et al. (2013, 2015), but contrasts with Theuerkauf and Rouys (2008); Figures 4,5,10

3. The distribution of red deer females was more strongly negatively correlated with wolves than that of red deer males, which is consistent with the selective killing of this sex by wolves in BF (Jędrzejewski et al., 2000; Jędrzejewski et al., 2002); Figure 4

4. The presence of spatial segregation between red deer males and females was in line with Kamler et al. (2008); Figure 8

5. The Strict Reserve was one of the hotspots of ungulate biomass (Jędrzejewska et al., 1997, 1994; Theuerkauf and Rouys, 2008); Figures 5,6,9,10

6. There were higher densities of red deer (especially of females) in the south part of the Strict Reserve, outside the wolf core area (Kuijper et al., 2015, 2013); Figures 5,6

7. There were higher detection rates of red deer and moose in canopy gaps (Churski et al., 2017; Kuijper et al., 2009); Figure 7

8. Wild boar had a preference for closed canopy deciduous dominated tree stands (Jędrzejewska and Jędrzejewski, 1998; Jędrzejewska et al., 1994; Kuijper et al., 2009); Figures 4,5,10

9. Ungulates had a general spatial association with deciduous (wild boar, bison) and mixed deciduous forests (red deer) (Jędrzejewska and Jędrzejewski, 1998; Jędrzejewska et al., 1994); Figures 4,5,10

10. Landscape use by bison was strongly influenced by supplementary winter feeding, the use of open areas inside the forest which are distributed sparsely within the study area and often associated with human settlements (Kowalczyk et al., 2013, 2011); Figures 4,5,6

All these, in addition to the formal methodological tests of our models, contribute to our confidence in the final results. We would like to repeat after (Grimm, 2005): “Useful patterns need not be striking; qualitative or’ weak’ patterns can be powerful in combination.”

4) Correlative evidence. It is important to acknowledge that your study draws inferences from correlative evidence, while establishing causality would require experimental manipulation of key parameters. Throughout the manuscript, terms are used that imply causality (to drive, determine, affect etc.), and these should be replaced with appropriate alternatives (to correlate, be associated with etc.) or combined with suitable qualifiers (may, seems, appears to be etc.). Two analyses in particular require more nuanced wording.

We agree about the general correlative nature of our results. We have carefully checked the manuscript and corrected the text where necessary, see our uploaded manuscript with tracked changes.

First, the 'herbiscape' analysis is identified as a primary point of novelty, and while conceptually this is correct, the evidence provided is largely correlative in nature. In the absence of direct experimental evidence of the effects of herbivores on the specific plant species examined (through controlled exclosure experiments), it is impossible to conclude that the spatial distribution of herbivores is causing the observed vegetation patterns. We recognise that such evidence is difficult to obtain, and we are not suggesting that additional work is required, but we think that this issue needs to be acknowledged (e.g., in subsection “Spatial heterogeneity in the landscape distribution of large herbivores” of the Discussion).

We agree on this point. Although the results of our study are in general correlative and based on non-manipulative observations of the system via remote sensing tools and field measurements, we have built our inference using other studies from Białowieża Forest, including our previous experimental studies on herbivory effects on vegetation (Churski et al., 2017; Kuijper et al., 2010) and quasi-experimental studies on the spatial patterns of ungulates browsing in relation to perceived risk effects (Kuijper et al., 2013). The studies by Kuijper et al., 2010, and Churski et al., 2017, were controlled exclosure experiments that provided the direct experimental evidence of the effects of herbivores on the specific plant species we examined in the current study (i.e. Acer platanoides and Carpinus betulus). The observed compositional changes in regenerating trees across herbiscapes were in accordance with these earlier experimental studies. Please, see the Discussion section for more details (subsection “Landscape scale herbivory regimes and their potential consequences for vegetation”). We believe that in combination with these previous exclosure studies our results allowed us to (at least) discuss the potential causalities as suggested by the observed patterns in the data and modelled correlations. However, at the same time we do not claim that the fine-scale patterns observed in the exclosure studies automatically imply corresponding causality at the landscape scale. Therefore we toned down our statements in the entire text, see uploaded manuscript with tracked changes.

Second, while it is clear from the data that red deer did not use the same areas as wolves, it is not possible to infer that they 'avoided' those areas (which implies causality). In fact, wolves may only select areas of high red deer density during periods of hunting; the negative correlation between deer and wolf distributions could therefore arise if areas used by wolves when they are not hunting are not preferred habitats of deer. Besides, as red deer are the major prey of wolves in the study area, why would wolves then not select areas of higher red deer density? We appreciate that experimental manipulation of this system is impossible, but note that analyses of spatio-temporal patterns (e.g., from simultaneous radio- or GPS-tracking of wolves and deer) would have considerably strengthened inferences. Please tone down your discussion of wolf-deer relationships accordingly.

We agree that our results are correlative in nature and do not present causality in the deer-wolf relationship. Thus, following the reviewer’s request we toned-downed our discussion accordingly. However, we believe, our camera-trapping system with random distribution and 24-h operation allows the monitoring of wolves both during hunting and non-hunting periods.

Answering the question of the reviewer “why would wolves then not select areas of higher red deer density?” we can cite her/his own answer: “In fact, wolves may only select areas of high red deer density during periods of hunting”. The reported average daily movement distance of wolves (based on telemetry) is 22-27 km (Jędrzejewski et al., 2001) and their average annual territories cover an average of 200 km2 (Jędrzejewski et al., 2007), which makes all our indicated red deer hotspots (based on camera trapping) within the daily range and within the territory of the present wolf packs. We would argue that wolves select breeding dens mainly based on human-related factors (in line with Sazatornil et al., 2016. Biological Conservation 201:103–110. DOI: 10.1016/j.biocon.2016.06.022) and then go to hunt in areas rich in deer.

The red deer is present across the entire landscape, although at much lower densities in the wolf high-use areas. In our opinion, it is very unlikely that this pattern is driven by the habitat preference of red deer, as suggested by the reviewer, as the wolf high-use areas also contain high-quality deciduous dominated and mixed deciduous forest stands intensively used by other ungulates (wild boar, bison, moose; Figure 2, Figure 5, also see Kwiatkowski, 1994). Covariates related to habitat quality were included in our models together with relative wolf encounter probability (our proxy for perceived predation risk). Our model results clearly point to the space use of the wolf being the main factor negatively influencing deer distribution rather than habitat factors. Whether these effects result from density- or behaviourally- mediated effects, or a combination of both, is less clear, as we discussed in subsection “The differential effect of large carnivores on the spatial structure of the large herbivore community”.

5) Presentation. This is a complex study, both in terms of the biological processes studied and the analysis techniques used, and it is critical to help readers follow the narrative. This is especially important as eLife's readership is very broad, and papers should be intelligible to non-specialists. To improve the presentation of the study, we suggest three main revisions. a) Technical terms should be avoided where possible, or at least clearly defined at the outset. For example, the Introduction talks about "cascading effects", "hyper-keystone species" etc. While these may be frequently-used terms in this research field, without further explanation, many non-specialist readers may struggle to understand the intended meaning. In general, the Introduction could be more focussed and reader-friendly.

We carefully edited the main text and added all missing definitions of more technical terms. We made sure that the main text is not overloaded with a technical jargon and added all missing definitions of more technical terms.

b) Since the Results and Discussion sections precede the Materials and methods, some aspects of the study system and methodology remain unclear until much later in the narrative. For example, the species under investigation haven't been mentioned until they are referred to in various parts of the Results section, and 'herbiscapes' are discussed in the Results before the concept is explained in the Discussion. There are two possibilities for dealing with this. Either the Materials and methods are moved forward, for a traditional I-M-R-D structure, or the front end of the narrative is carefully revised, to equip the reader with all the information they need to understand the results and their interpretation.

We moved the methods forward as suggested.

c) We think it would really help readers if there was a 'graphical abstract' that illustrates, in two separate panels: (i) how data were collected and processed; and (ii) the biological relationships uncovered by the analyses (humans --> wolves --> prey --> vegetation). Panel (ii) could provide explicit call-outs to figures and materials in the online supplement that provide support for a given link, helping readers navigate these various materials.

We prepared the’ graphical abstract’ as requested. We want to thank the reviewers for this suggestion! We think that it was a great idea and we had a lot of fun working on it together with a specialist in this field and our friend Lisa Sanchez.

6) Human activity. Did your camera traps record human activity (on/off established tracks)? If so, why were these data not included in the models, instead of (or in addition to) your chosen proxy variable (distance to settlements)? Appendix 1—figure 3 shows human activity data, but it is unclear how these were obtained, and "combined" with other datasets, and whether any of this was used in the models. Incidentally, human detections should enable a very useful assessment of the reliability of your camera-trapping method (the distribution of human observations is expected to match the known distribution of human activity in the forest); see also comment #3 above.

Yes, our camera traps recorded humans and we originally intended to use this data as a variable representing the human impact. However, it appeared that the design of both our camera trapping surveys (for ungulates and carnivores) precluded obtaining an unbiased measure of human activity over the BF landscape. A preliminary exploration of the raw data showed a very random pattern of the human record distribution in both surveys, and we considered it unlikely to reliably reflect the true spatial pattern of human interference in BF. In the ungulate survey, the cameras were placed pseudo-randomly in the forest and within a distance of 100 m from the nearest roads (both paved and unpaved), large clearings and settlements, and provided very few human records. In contrast, the carnivore survey, with cameras placed on forest roads, provided more abundant human data, but with records biased heavily by vehicle traffic. The reason for this lies partly in our selection of forest roads as locations for camera traps – we used both very small forest roads not utilized by cars and larger roads used e.g. by forestry transportation machines and border guards, who patrol the area regularly in their vehicles. Also, we did not have cameras monitoring the major (paved) roads in BF that are open for the public and cars (see Figure 1) (i.e. the roads with heavy traffic in the study area).

For these reasons, we decided to express the landscape scale human impact using spatial covariates derived directly from the available GIS data. We believe that the chosen distance-based covariates (distance to all settlements and distance to major settlements) represent human (especially non-motorised) pressure well, as they consider permanent human-made structures and the resulting enduring effect of human interference with wildlife. The distance to major roads was highly correlated with the distance to major settlements (0.87; Appendix 1—figure 1). Another proxy describing the landscape-scale variation in human pressure was the density of protected areas, which are areas of minimal human activity due to exclusion of most forestry practices, hunting and motorized vehicles. We believe that the combination of these variables provided more accurate information on landscape scale and long term human pressure on wildlife populations than our camera trapping surveys, which were not specifically designed for the purpose of collecting this kind of data. As a visualization we included the available data from STRAVA (the global “heatmap” of human activity, https://www.strava.com/heatmap#10.75/23.64564/52.74927/hot/all), see Appendix 1—figure 3. The obtained picture of human activity (people using the STRAVA mobile application) fits well to our covariates describing landscape scale human impact.

7) Sampling periods. The analyses assessing patterns for a subsampled ungulate dataset (for the months covered by the carnivore survey) are very useful. But, it is not clear why an extra month was included here (Aug-Oct; subsection “Carnivores survey”). Why not use Sept-Oct, to make it strictly comparable? More importantly, you say that the "results were similar to those based on the full dataset", yet careful comparison between Figures 4/5 and Appendix 1—figures 36/37 reveals important differences. This should be acknowledged, further investigated, and discussed. Which of your main conclusions remain unchanged when you use this subsampled dataset, and which ones are affected?

The differences between the results based on the full dataset and its subset covering only three months was to be expected. With the full dataset covering all seasons we obtained a picture of the yearly-averaged responses of ungulates to factors affecting their distribution over the landscape. It is to be expected that responses can somewhat differ for some species and factors when we limit the observation window to a short period of only 3 months. Considering the logistics of the entire field work we planned to assure a good, close-to-uniform spatial coverage of the entire study area for the entire study period (2 years). Additionally, we aimed to ensure a similar level of spatial effort for two distinct seasons, i.e. green-on and green-off. Thus, a three-month-only subset of data would not have been representative in terms of the spatial coverage. From a technical point of view, our spatial models are (very) data hungry and we thus chose to use the full dataset, which provided a larger sample size and better spatial coverage of the study area, both of which are needed for making robust inferences using complex hierarchical spatial models such as ours (see Statistical model section). For these technical reasons we included three months in this extra analysis to ensure convergence of the models. When limiting ourselves to only two months, we had too little data to fit the models, which did not converge.

In our opinion, this extra analysis is not critical for confidence in our analysis and results (see our responses to points #1 and #3), especially when taking into account the fact that our model of wolf space use matches well the telemetry and den location data from the past (see our response to point #1). All this indicates the high level of stability in the way that wolves use the Białowieża Forest landscape.

8) Broader implications. The study makes a notable contribution to our understanding of key processes in community ecology, but we felt it would be useful to discuss some broader implications in the Conclusions section. Specifically, are any of the findings relevant for ongoing conservation/management efforts, in the study area or elsewhere?

Following this suggestion we have added a new paragraph to the Discussion. Please see paragraph six of subsection “Landscape scale herbivory regimes and their potential consequences for vegetation”.

9) Transparent reporting. Some of the points raised in eLife's guidance document require attention (e.g., power analyses, which relates to comment #1 above).

Please see our response to point #1. Also, we are clear and transparent about this issue in the submitted transparent reporting form:

“No explicit power analysis was used prior to our field data collection using camera traps. As our study was of an observational, rather than experimental nature, and as we investigated complex spatial interactions between many species representing multiple trophic levels, we were unable (to our knowledge) to run an explicit power analysis. However, in our methods we ensured a large enough number of sampled locations to successfully fit our complex spatial models (diagnostics are presented in Appendix 1). Our sampling methodology, together with a short discussion of all the arbitrary decisions we had to take, are described in detail in the main text in the Materials and methods section.”

We would also like to refer to the work by Parsons et al., 2017, recently published in eLife: “No explicit power analysis was used to predetermine sample size. Our sample size goal was 20 spatial replicates (equating to ~420 trap nights), which has been found to maximize precision for estimating detection rate (Kays et al., 2010; Rowcliffe et al., 2008). Camera sites are biological replicates, parallel measurements capturing random biological variation. This study did not include technical replicates.”

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The reviewers have carefully evaluated your revised manuscript and response letter. While reviewer #2 was happy with the revisions overall, apart from a few minor points that still need addressing, reviewer #3 remained concerned about data quality, and specifically, the mismatch in timing of the predator and prey datasets. The reviewers' original reports are appended below for reference. We share reviewer #3's concerns, but are prepared to consider a final revision that addresses these points more explicitly and tones down claims further. It would also be important to better acknowledge throughout that correlations between spatial species distributions were examined, not interactions (as for example implied by the use of the term “cascading effects”).

Reviewer #2:

All of my major and minor comments have been addressed in this revised version of the manuscript. I believe the revision is much-improved and would recommend publication after addressing the few minor comments below.

Minor Comments:

Introduction paragraph two: Unsure of the accuracy of this claim that this subject has often been ignored – "The strong spatial variation in herbivory impact, resulting from the space use of by different functional groups of herbivores has often been ignored." Suggest re-wording to simply state the effect, as it has been shown elsewhere, e.g. "There is often strong spatial variation in herbivory impact resulting from the space use of by different functional groups of herbivores."

We have changed this sentence according to the reviewer's suggestion.

Subsection “Spatial heterogeneity in the landscape distribution of large herbivores” paragraph three: Avoid repeated and strong claims of novelty (as this type of study has been done in other systems, as you say). Suggest re-wording to "Similar studies in African ecosystems have shown that…"

We have changed this sentence according to the reviewer's suggestion.

In the same subsection: Suggest rewording to focus specifically on the benefits of this study from providing information on a temperate ecosystem – this is the major point of novelty and, in my view, what makes this study a valuable contribution to this body of literature, and should be emphasized. The current wording does not make this clear.

We have reworded this sentence according to the reviewer's suggestion:

“Our contribution to all the above studies is a high-resolution picture of the spatial structure of an entire community of large herbivores that incorporates both bottom-up and predation- and human-related top-down factors. We believe one of our major points of novelty is in providing information on these spatial interactions for a temperate ecosystem.”

Reviewer #3:

The first concern regarding the manuscript was the robustness of statistical inference. In the rebuttal, the authors assume that if there is no problem with the robustness of the statistical analyses, then the results are reliable. However, the problem does not lie with the analyses but with the quality of data. If the data are not representative in the first place, then there is even no way to verify the robustness of the results. The major problem I see is the temporal mismatch of data sampling. While you sampled ungulate distribution over 2 years, you only sampled predators during a single period of 2 months. Besides, you state that "In August-September, pups begin to travel with other pack members and move more widely through their territory, returning to the core of their territory on a regular basis (Jędrzejewski et al., 2001)". This is however not completely correct. Wolves actually use core areas only during the pup rearing season but not during the rest of the year. Therefore, seasonal data is not representative for the year-round space use. You can therefore compare space use by wolves in September-October only with prey space use during the same season. Otherwise, it is possible you will get random results. Even if prey space use in September-October is representative for the whole year (using the whole year will actually create only the impression of more robust data by artificially increasing sample size, but as the predator sampling is reduced to a short period, it does not help), wolf space use in September-October is definitely not representative for the whole year and therefore also not the best period to quantify space use as you state at the end of the first paragraph in subsection “Carnivores survey”. I would therefore not agree with your statement in the rebuttal regarding seasonal differences in wolf space use: "We appreciate this idea but we did not have data on large carnivores from other seasons. Moreover, we expect little variation in wolf use of the landscape". On the contrary, there are huge differences in wolf space use comparing spring/summer and autumn/winter. This is especially important as September-October also overlap with the rut of red deer, thus the period when wolves rather select hunting male and not female deer (during this period you found that females more correlated with wolves). Given that you did also not consider any temporal information (times of hunting vs times of resting), it is not really possible to assess if the correlations are actually based on real predator-prey interactions or if they are random results. It is possible that the results actually are results of trophic relations, but the data structure does not allow explicit conclusions.

The critique of reviewer #3 essentially comes down to the question of whether the wolf data collected by our camera-traps in September-October is representative of the year-round and landscape-scale space use of wolves. The reviewer states that it is not, and that the temporal mismatch of camera trap-based sampling of large carnivores (in autumn) and ungulates (year-round) could ultimately lead to “random results”. The reviewer specifically argues that “wolves actually use core areas only during the pup rearing season but not during the rest of the year”. Below we argue that our assumption is defendable and that our camera-trapping carnivore-data can be regarded representative:

1) Firstly, the reviewer’s statement that wolves use core areas only during the pup rearing season but not during the rest of the year is contrary to extensive earlier studies from the study area. The year-round, rotational use of core areas by wolves, including outside the pup-rearing period, was shown by Jędrzejewski et al., 2001, on the basis of extensive radio-telemetry work (Figures 2 and 3). The presence of clearly defined core areas during bimonthly intervals throughout the year was also shown by Jędrzejewski et al., 2004, in their detailed study on the process of wolf pack splitting. Both Jędrzejewski et al., 2007, and Schmidt et al., 2009, showed that on a yearly basis core areas are the territory areas used most intensively by wolves. Furthermore, Zub et al., 2003 also showed that in winter (study period from November till April) the density of territorial markings is concentrated in core areas due to more intensive use of core areas by packs. Finally, the study of Nowak et al., 2007, found that spontaneous howling activity in summer and late autumn is strictly connected to the core areas of their territories, where the breeding den had been located, and not from the peripheries. In summary, the existing literature clearly shows that core areas are the most intensively used parts of the annual wolf territory throughout the year and many wolf activities additional to breeding and pup-rearing concentrate in this area.

That wolves continue to use the core area of their annual wolf territory also after pup rearing, although with different intensity, was also shown by Theuerkauf et al., 2003. The core areas are highly selected in both summer and winter (see Figure 2 in Theuerkauf et al., 2003). We generally agree with the reviewer that there are seasonal changes whereby the summer season (including the pup-rearing period) is characterized by a stronger concentration in the core areas. However, these seasonal changes in territory use do not change the overall pattern of concentrated use in some parts of their territories.

2) The seasonal changes in wolf movements are most clearly marked in adult breeding and sub-adult females. Indeed, during the pup-rearing period their movements are restricted to the vicinity of the den. However, other members of the pack (adult males, adult non-breeding females, subadult males) do not show such “huge differences”. Please also note, that camera traps collect data on wolf presence and activity randomly, irrespective of the sex or breeding status. The period from September to October was chosen to monitor the activity of carnivores, considering that both wolves (of all reproductive status) and lynx then resume travelling long distances within their home ranges after the breeding season. According to previous telemetry studies conducted on both predators in the study area wolves covered the largest distances and ranges during this period as compared to other months (Jędrzejewski et al., 2001), and lynx (all except breeding females) were equally mobile year-round (Jędrzejewski et al., 2002). We thus believe that as we were unable to conduct a year-round survey of carnivores, September – October is the most suitable period for sampling the activity of carnivores at the scale of the whole forest complex as it reflects the general annual activity of large carnivores most reliably.

3) In our previous response to this particular critique, we already showed that the non-uniform space use of wolves as predicted by our camera trap study, is in close correspondence with earlier radio-telemetry studies (showing the year-round space use of wolves) and with the distribution of known locations of breeding dens (see Appendix 1—figure 4). These independent data sources all show that the landscape use of wolves is neither uniform nor random and core areas (with a higher intensity of use) are clearly “visible” at the landscape scale. Moreover, the core areas and spatial arrangements of territories of the four wolf packs in Białowieża forest (BF) have been relatively stable over the last decades. We elaborated on this in the main text (subsection “Specification - ungulates model” in the Materials and methods and subsection “Landscape scale herbivory regimes and their potential consequences for vegetation” in the Discussion) and in our previous response letter (points #1c).

4) We showed in our previous response letter that restricting our ungulate dataset to only part of the year, to create a better temporal match between wolf and ungulate data, does not change the model outputs and conclusions of our study related to large carnivores (see Appendix 1—figures 36 and 37).

Conclusions:

We do not disagree with the reviewer that wolves show seasonal changes in their territory use, but the arguments provided above illustrate that there is no rationale for expecting these changes to undermine the patterns we found in the study. Based on existing knowledge and our data, we can expect continuous higher relative use of core areas than of other parts of their territories. Thus the spatial patterns of wolf landscape use observed in autumn, when they clearly concentrate in their core areas but also widely use their entire territories, gives a representative, ‘averaged’ picture of the relative use of the BF landscape throughout the year. Therefore, we feel confident that our collected wolf data is representative for showing year-round patterns of wolf space use, despite the seasonal variation in territory use over the year (see for example Theuerkauf et al., 2003) referred to by the reviewer.

Our analyses show that these core areas have an effect on the space use of red deer. This does not exclude the possibility that additional finer-scale and more direct interactions occur between wolves and their prey during instances of actual wolf hunting activity. Much of the hunting activity takes place outside the wolf core areas and all red deer hot spots predicted by our models are within the daily range of each wolf pack in BF. We elaborate on this topic in our previous response to the reviewers (see point #4b). However, our data suggest that the continuous presence and higher relative use of wolf core areas has a landscape-level effect on the distribution of their main prey species (red deer) on a yearly basis.

If the reviewer continues to disagree with the arguments we have provided, we feel it may be more appropriate to continue this discussion on the basis of a scientific comment to our study that we would be happy to reply to.

References:

Jędrzejewski W, Schmidt K, Theuerkauf J, Jędrzejewska B, Kowalczyk R. 2007. Territory size of wolves Canis lupus: Linking local (Bia?Owie?A Primeval Forest, Poland) and Holarctic-scale patterns. Ecography30:66–76. doi:10.1111/j.2006.0906-7590.04826.x

Nowak S, Jędrzejewski W, Schmidt K, Theuerkauf J, Mysłajek RW, Jędrzejewska B. 2007. Howling activity of free-ranging wolves (canis lupus) in the białowieża primeval forest and the western beskidy mountains (poland). Journal of Ethology25:231–237. doi:10.1007/s10164-006-0015-y

Zub K, Theuerkauf J, Jędrzejewski W, Jędrzejewska B, Schmidt K, Kowalczyk R. 2003. Wolf pack territory marking in the Bialowieza Primeval Forest (Poland). Behaviour140:635–648. doi:10.1163/156853903322149478

https://doi.org/10.7554/eLife.44937.062

Article and author information

Author details

  1. Jakub Witold Bubnicki

    Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    kbubnicki@ibs.bialowieza.pl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2064-3113
  2. Marcin Churski

    Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland
    Contribution
    Conceptualization, Supervision, Validation, Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8727-0203
  3. Krzysztof Schmidt

    Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland
    Contribution
    Conceptualization, Resources, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9043-291X
  4. Tom A Diserens

    Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland
    Contribution
    Writing—review and editing
    Competing interests
    No competing interests declared
  5. Dries PJ Kuijper

    Mammal Research Institute, Polish Academy of Sciences, Białowieża, Poland
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared

Funding

National Science Centre, Poland (2012/07/N/NZ8/02651)

  • Jakub Witold Bubnicki

National Science Centre, Poland (2015/17/B/NZ8/02403)

  • Marcin Churski
  • Dries PJ Kuijper

EuroNatur (PL-15-500-28)

  • Krzysztof Schmidt

National Science Centre, Poland (2017/25/B/NZ8/02466)

  • Tom A Diserens
  • Dries PJ Kuijper

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Roman Kozak, Andrzej Waszkiewicz and Tomasz Kamiński for their field and lab work, and to Lisa Sanchez for preparing a graphical abstract of the article. We thank Bogumiła Jędrzejewska for her contribution at the initial stage of the project. We thank the administration of the Białowieża National Park, the administrations of Białowieża, Hajnówka and Browsk Forest Districts, and Polish authorities (Ministry of Environment, Polish General Directorate of Environment Conservation and Regional Directorate of Environment Conservation in Białystok) for permissions to work in Białowieża Forest. We thank European Space Agency (ESA) and Planet Labs Germany GmbH for providing RapidEye satellite imageries (ESA TPM PI project no. 28377). LIDAR and vegetation data used in the analysis were kindly provided by the team of the ‘LIFE+ ForBioSensing’ project which is co-funded by the EU Life Plus programme (contract number LIFE13 ENV/PL/000048) and the National Fund for Environmental Protection and Water Management in Poland (contract number 485/2014/WN10/OP-NM-LF/D). The work of JWB was supported by the National Science Centre, Poland (grant no. 2012/07/N/NZ8/02651). The work of DPJK and MC was supported by the National Science Centre, Poland (grant no. 2015/17/B/NZ8/02403), and the work of KS by European Nature Heritage Fund EURONATUR, Germany (grant no. PL-15-500-28). In addition, the work of DPJK and TAD was supported by funding from the National Science Centre, Poland (grant no. 2015/17/B/NZ8/02403 and grant no. 2017/25/B/NZ8/02466.

Senior Editor

  1. Ian T Baldwin, Max Planck Institute for Chemical Ecology, Germany

Reviewing Editor

  1. Christian Rutz, University of St Andrews, United Kingdom

Reviewers

  1. Christian Rutz, University of St Andrews, United Kingdom
  2. Joern Theuerkauf

Publication history

  1. Received: January 7, 2019
  2. Accepted: September 13, 2019
  3. Accepted Manuscript published: October 2, 2019 (version 1)
  4. Version of Record published: October 22, 2019 (version 2)
  5. Version of Record updated: November 19, 2019 (version 3)

Copyright

© 2019, Bubnicki et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 833
    Page views
  • 136
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)