Highlights

  • Ungulate presence cascades their effect to lower trophic levels.

  • Decomposition and mineralization rate remain unchanged due to ungulate presence.

  • Ungulate abundance strengthens the cascading effects in this temperate system.

  • Ungulates may trigger slow feedback loops on nutrient cycling and forest succession.

Introduction

Through the input of nutrients via urine and feces, as well as trampling and selective browsing, ungulates can modify above-ground vegetation structure and composition (Ramirez and others 2019), which in turn can trigger below-ground large cascading effects, including soil quality, invertebrate composition and decomposition rates (Allombert and others 2005b; Bressette and others 2012). This may trigger a feedback loop, in which the altered nutrient availability affects the abundance and composition of trees and ungulates (Wardle and others 2002). Ungulate feeding preferences theoretically determine the speed that the feedback loop acts: Preferential feeding on slow-growing species with nutrient-poor leaves results in a fast feedback loop and a more fertile and productive ecosystem, whereas preferential feeding on fast-growing species with nutrient-rich leaves results in a slow feedback loop and a more infertile and unproductive ecosystem (Ritchie and others 1998; Wardle and others 2002; Andriuzzi and Wall 2017). Yet, empirical evidence suggests that herbivory can modify the amount of carbon allocated below-ground by plants (Kardol and others 2014), which also affects the speed of the feedback loop.

Feedback loops may become more important with the loss of apex predators in temperate forests (for example, wolves) because ungulates are released from top-down control, resulting in higher ungulate abundance. Hence, higher numbers of ungulates increase pressure on lower trophic levels and the potential for trophic cascades and feedback loops (Kardol and others 2014). However, important literature advocates that feedback loops across biomes are strongly mediated by primary productivity, seasonality, natural disturbance (for example, fires and storms) and the ratio between conifers and broadleaves (Pastor and others 1988; Hobbs 1996; Augustine and McNaughton 1998; Frank and others 2000; Pringle and others 2007). Thus, feedback loops are complex mechanisms that are highly dependent on multiple external factors.

In any ecosystem, ungulates may indirectly affect decomposition via changing litter quantity and quality, which affects the abundance of detritivores (Moretto and others 2001; Mason and others 2010; Frouz and others 2015). Ungulates can also reduce litter depth on the forest floor by mixing the litter with soil as they scrub the ground with their hoofs (Hobbs 1996) and by browsing on broadleaf species, which subsequently changes the species composition of the stand to evergreen coniferous species that have a lower litter production rate (Pastor and others 1993; Husheer and others 2005; Ramirez and others 2018). Alternatively, ungulates may increase litter depth by shifting species composition toward conifers which have more recalcitrant litter, resulting in an accumulation of litter over time. The relative importance of these two processes determines in the end the depth of the litter layer.

Invertebrates can regulate nutrient cycling and soil structure by feeding, nesting and burrowing on leaf litter (Jones and others 1994; Brussaard 1997; Lavelle and others 2006; Kamau and others 2017). Detritivores decompose litter first by fragmenting it into small pieces so bacteria, fungi and protists can then convert organic compounds into inorganic compounds like phosphate, ammonium, water and carbon dioxide (Aerts 2006). Decomposition is also highly dependent on the micro-environment (Vasconcelos and Laurance 2005). For instance, ungulates can increase soil compaction and light availability in the understory (Ramirez and others 2019), thus triggering cascading effects to other trophic levels. Soil compaction leads to a reduction in soil porosity, which in turn results in less oxygen and resources that are available for invertebrates, which ultimately slows decomposition rates (Lavelle and others 1992; Van Klink and others 2015). Overall, trampling may reduce the amount of soil biota (Bressette and others 2012) and thus slow down decomposition by invertebrates, fungi and bacteria (Hättenschwiler and others 2005).

Several studies have evaluated cascading effects of wild ungulates on temperate forests in Europe, North America and Oceania (Wardle and others 2002; Allombert and others 2005b; Bressette and others 2012; Kardol and others 2014; Lilleeng and others 2018) by using fencing experiments that included or excluded ungulates. These studies suggest that there are strong top-down effects of ungulates on vegetation composition, which in turn influence soil invertebrate communities and decomposition rates (Allombert and others 2005b). However, none of these studies have evaluated the cascading chain across understory vegetation, soil attributes, soil invertebrates, litter decomposition and nutrient mineralization. For example, most studies only measured the effects of herbivores on invertebrate composition (Allombert and others 2005b; Bressette and others 2012; Lilleeng and others 2018) or decomposition rate (Wardle and others 2002; Andriuzzi and Wall 2017). Hence, there is a need to evaluate ungulate cascading effects across the complete chain.

To understand the direct and indirect effects of ungulates on temperate forests, we conducted a study in which red deer (Cervus elaphus), roe deer (Capreolus capreolus), fallow deer (Dama dama), wild boar (Sus scrofa) and mouflon (Ovis orientalis) were experimentally excluded. We compared vegetation, litter layer, soil quality, rodent activity, leaf litter invertebrates, leaf litter decomposition and nutrient mineralization between paired fenced and unfenced plots established in twelve forests across the Veluwe, the Netherlands. We evaluated litter decomposition using litter of three tree species that dominate forest succession on sandy soils (Betula pendula, Pinus sylvestris and Quercus robur) and nutrient mineralization with resin bags. We predicted that ungulates, by means of browsing and trampling, can trigger cascading effects on lower trophic levels, including litter decomposition and mineralization rates (Figure 1), by restricting tree diversity and density, litter depth, invertebrate diversity and biomass as well as soil quality (McGarvey and others 2013; Lilleeng and others 2018; Ramirez and others 2019). Due to feedback loops, we theoretically expect a decrease in tree diversity that regenerates in unfenced compared to the fenced plots because ungulate browsing and trampling reduces nutrient cycling in soil, thus enhancing bottom-up control of vegetation (Ritchie and others 1998; Wardle and others 2002).

Figure 1
figure 1

Conceptual model of the cascading effects (represented by solid lines) of wild ungulate presence, treatment age and ungulate abundance on forest regeneration (diversity and density), litter (depth and quality) soil (compaction and pH), rodent activity, invertebrates (diversity and biomass), litter decomposition (tea, pine, oak and birch) and nutrient (nitrogen and phosphate) mineralization in soil. Feedback loops (not tested in this study) are represented with dashed lines. The relationship between components are represented by lines, and the symbols positive (+) and negative (−) indicate the direction of the relationships.

Full size image

Methods

Study Area

Data were collected at the Veluwe, a 1200 km2 area located in the central part of the Netherlands (52° 5′ N, 5° 48 ′E). Mean annual precipitation is 900 mm y−1 and mean annual temperature is 9.4°C (Kuiters and Slim 2002). The main soil types are xeric humic podzols and brown earths (inceptisols), depending on the parent material that ranges from aeolic driftsand and cover sands to Pleistocene loamy fluvioglacial sands (Kuiters and Slim 2002). The Veluwe is covered by a mosaic of forests, drift sands and heathland, where forests cover two-thirds of the total area. The main tree species are Scots pine (Pinus sylvestris), European oak (Quercus robur), European beech (Fagus sylvatica), Japanese larch (Larix kaempferi), Douglas fir (Pseudotsuga menziesii) and European birch (Betula pendula).

Historical records of ungulate abundance and composition at these sites are not available; hence, we obtained current ungulate abundances from unpublished material. This was done by deploying 21 camera traps mounted to trees at 50 cm height, with a 70 m spacing, within 1 km2 of the forest surrounding each of the fenced plots. The deployment took place in the summer of 2017 with a total of 450 camera days−1. As a result, we were able to quantify ungulate trap rate for each of the experimental sites. The main species recorded were roe deer, fallow deer, red deer and wild boar, with an average density of 13.6 animals per km−2 in 1998 (Kuiters and Slim 2002), and the composition varied across sites.

Experimental Design

Research was carried out at twelve sites distributed across the entire Veluwe region (Figure 2, Online Appendix A.1), in an area of 70 by 40 km. At each site, a single fenced forest plot had been established to protect forest regeneration from ungulate browsing and was paired with an unfenced control plot, approximately 10 m apart. For plots without a control plot, we randomly assigned a control area within 10 meters distance from the fenced plot. Plot size was on average 631 m2 (range 100–2100 m2), and fencing consisted of a 2.10 m tall metal fence. The fenced/unfenced plots were established 1–33 years ago, in logging gaps. Because historical data are not available on gap size and the fact that forestry in the Netherlands promotes forest regeneration with small gaps (30–50 m in range), we assume that all gaps were the same size.

Figure 2
figure 2

Map of the twelve research sites at the Veluwe, the Netherlands. The dark gray area is covered by forest. 1. Hoenderloo, 2. Gortel, 3. Ullerberg south, 4. Ullerberg north, 5. Hoge Veluwe north, 6. Hoge Veluwe south, 7. Achterpark, 8. Garderen, 9. Dellen, 10. Epe, 11. Oostereng and 12. Rheden (Online Appendix A.1).

Full size image

Above-ground Vegetation

To quantify vegetation structure and diversity, the vegetation was surveyed during spring and summer of 2016 and 2017. Within each pair of fenced/unfenced plots, 5 × 5 m quadrats were randomly established by drawing numbers for the x and y axis, which represented a coordinate system. We established two quadrats per plot when regeneration heterogeneity was low (that is, low species diversity and little variation in forest structure, which is typical for smaller fenced/unfenced plots), and three or four quadrats per plot when heterogeneity was high. Data from all 2–4 quadrats were averaged to obtain a single value per plot. Each woody individual above 15 cm within the quadrat was identified to species level and measured for its height.

Relative abundance and Shannon diversity were then calculated for all tree saplings (Shannon 1948). Stem density was quantified as the number of stems per quadrat and then down scaled to m2. Ungulate browsing was scored as “0” if no shoots were damaged by browsing or “1” if one or more shoots were browsed. Plot litter quality was derived by scoring all saplings (> 15 cm in height) with the species litter quality score as in Maes and others (2019) and then averaging the values for the entire plot. Litter quality scores range from 1 to 5, where high-quality litter scores 5 and low-quality scores 1. This was done to quantify the quality of litter available for decomposition. We measured litter depth, because litter can act as a barrier for the establishment of woody species from the seedbank (Facelli and Pickett 1991; Schramm and Ehrenfeld 2010), while it can provide shelter and resources for invertebrates. Litter depth was measured with a ruler at two points in each quadrat, as the distance between the litter surface layer and the mineral soil. Finally, rodents are important seed predators and dispersal agents in temperate forests (Nathan and Muller-Landau 2000; Iida 2006). Thus, rodent activity was quantified as the number of litter and tea bags that were removed (Mdangi and others 2013); this assumption was confirmed by field observation.

Decomposition and Mineralization Rate

To quantify litter decomposition rates, we used litter bags with natural leaf litter of three dominant tree species and tea bags as a standardized litter. The three species were the coniferous Pinus sylvestris (pine), and the broad-leaved Betula pendula (birch) and Quercus robur (oak). These species were selected because they occupy different successional positions in the forest and their leaves have different nutrient content which makes them more or less palatable for ungulates (Maes and others 2019). In each plot, 18 leaf litter bags were buried on the ground, belonging to the three species (6 bags per species). Litter bags consisted of 15 × 15 cm polyethylene mesh with a pore size of 1.5 mm, filled with approximately 2.0 grams of air-dried litter (Gartner and Cardon 2004; Bärlocher 2005). Additionally, twelve tea bags were placed in each plot: six with Lipton Rooibos tea and six with Lipton Green tea (Keuskamp and others 2013), in order to have a decomposition base line with standardized litter. Furthermore, by using two tea types with contrasting decomposability: green tea (fast decomposition) and rooibos (slow decomposition), a decomposition curve can be drawn from using a single measurement in time (Keuskamp and others 2013). To quantify nutrient mineralization, two resin bags per plot were prepared as in Lajtha (1988): with four grams (dry weight) of Dowex 50W-X8 cation-exchange resin and then buried below the litter layer to quantify nitrogen (N) and phosphate (PO4) mineralization rate in soil. During spring of 2017, when decomposition started, all tea, resin and litter bags were buried below the litter layer in order to prevent uprooting by ungulates. The tea bags and resin bags were harvested at 3 months because at this time their mineralization curve stabilizes (Keuskamp and others 2013), half of the litter bags were harvested after 6 months, and the remaining litter bags were harvested after 12 months in order to have decomposition constants at two moments in time. Harvested bags were oven-dried at 60°C for 3 days, dirt particles were gently removed with an air blow gun, and litter dry mass was measured. Resin bags were prepared and processed as in Lajtha (1988), and samples were analyzed with the segment flow analyzer.

Then, decay rates for leaf litter mass and soil nutrient mineralization were calculated using the negative exponential decay model (Hasanuzzaman and Hossain 2014):

$$ \frac{{\text{Bf}}}{{\text{Bi}}} = {\text{e}}^{ - k*t} $$

where (Bf) is the final biomass and (Bi) is the initial biomass. (k) is the decomposition rate and (t) is the time. The annual decomposition exponent (k) for pine, oak and birch litter was calculated as:

$$ k = 365*\frac{{\ln \left( {\text{Lf}} \right) - \ln \left( {\text{Li}} \right)}}{\text{time}} $$

Time was the duration in days that samples were in the field. Nitrogen and phosphate mineralization rates are calculated as:

$$ {\text{Nutrient}} \left( {\frac{{\text{mg}}}{{\text{l}}}} \right) = \left( {\frac{C}{{\frac{1000}{V}}}*\frac{1000}{W}} \right)/{\text{time}} $$

where (C) is the nutrient concentration in the soil extract, (V) is the volume of potassium chloride used for the extraction, (W) is the exact weight of the resin sample used for extraction, and (time) is the duration in days that samples were on the field (Griffin and others 1995).

Below-ground

To quantify ungulate trampling, we measured soil compaction. Soil compaction is defined as the densification of soil grains in a given area as a result of external pressure and it is important for nutrient cycling in soil (Bassett and others 2005). For each location where a litter bag was placed, soil compaction (in kPa) was measured with a soil penetrometer (Eijkelkampm Serie 2518611) (Bassett and others 2005). Because soil pH affects tree development by controlling plant nutrient availability (Lucas and Davis 1961), pH was measured by extracting a soil sample at a depth of 10 cm from each litter bag placement and measuring it with a pH meter (inoLab pH Plus). We quantified invertebrate diversity and biomass because invertebrates play an important role in litter decomposition. We extracted two soil samples of 2.00 kg from each plot, and with the use of a metal screen and tweezers, macro-invertebrates (> 1.5 mm) were collected and stored in 70% ethanol. All invertebrates were identified to order level with the use of a compound microscope (Novex, P-series) and a soil invertebrate guide (Krogh 2010). Invertebrate dry biomass volume was obtained by measuring with a digital caliper (Mitutoyo 500), with a resolution of 0.02 mm, the invertebrate length × width × height and controlled by the soil sample weight. With this information, invertebrate Shannon diversity was then calculated.

Data Analysis

All statistical analysis was done in R version 3.4.0 (Team 2013). To assess how the 15 response variables (Table 1) were associated among themselves and with fencing, treatment age and ungulate abundance, we carried out a principal component analysis (PCA) using the “vegan, version 2.5-2” package (Oksanen and others 2013). For this we used average data. To quantify the effect of fencing on each of the 15 forest response variables (individual data), we used general linear mixed models (GLMM) because all variables were continuous and normally distributed, with forest site as random grouping factor. Akaike’s information criterion (AIC) was used to select the best-fitting model. The best fit is all models within 2 AIC units from the lowest AIC. This was done by using the packages “nlme, version 3.1-137” and “stats, version 3.4.0” (Pinheiro and others 2014).

Table 1 General Linear Mixed Model (GLMM) Results for Treatment Effect (fencing = 1 vs. unfencing = 0) on 15 Response Variables
Full size table

To test our conceptual model of ungulate cascading effects on lower trophic levels (Figure 1), a pathway analyses model (PAM) was used (Shipley 2016; van der Sande and others 2017). PAM is a multivariate statistical method that tests direct dependencies among a set of variables in complex path networks and normally consists of dependent and independent variables. A total of 15 models were needed to test all the pathways of the conceptual model, and each model consisted of one response variable and one to five predictors, depending on the tested path section. For example, litter depth depends on ungulate presence, sapling density and sapling diversity, whereas soil compaction only depends on ungulates. Forest site was included as random factor to statistically control for variation across locations. All data used in the models were averaged. The packages “lavaan, version 0.6-3” and “lavaan.survey, version 1.1.3.1” were used for this analysis (Rosseel 2012). Since regeneration, lower trophic levels do not only respond to ungulate presence, we conducted GLMM’s with treatment age and ungulate abundance as predictors and regeneration, litter and soil as response variables (averaged data), with forest site as a random grouping factor. We also tested for the interaction between ungulates and treatment and the results indicated non-significance; hence, we opted not to include it in the analysis. AIC values were used to select the best-fitting model.

Results

Treatment Effect: Fenced Versus Unfenced

The PCA analysis indicated clear correlations among the 15 response variables (Figure 3). Principal component 1 (Dim. 1—horizontal axis) explained 25.7% of the variation and was mainly associated with sapling density, sapling diversity, rodent activity, litter depth, tea decomposition, invertebrate diversity and invertebrate biomass at the right and soil pH and soil compaction at the left side. Principal component 2 (Dim. 2—vertical axis) explained 18.9% of the variation and was associated with litter quality, oak, birch and pine decomposition at the upper section and with nitrogen and phosphate mineralization at the lower section. Ungulate fence effect was mainly found along the second PCA dimension: Fenced plots without ungulates were found at the right side (with the 95% confidence interval indicated by a red ellipse) and unfenced plots with ungulates were found at the left side (indicated by a green ellipse). The ellipse of the unfenced treatment was much smaller than the ellipse of the fenced treatment, indicating that ungulates had a homogenizing effect on lower trophic levels. The twelve forest sites separated more out on the second dimension, indicating that site conditions modify lower trophic levels.

Figure 3
figure 3

Principal component analysis (PCA) for 17 variables related to treatment, treatment age, ungulate abundance, vegetation, soils, decomposers, decomposition and mineralization rates in forests across the Veluwe, the Netherlands. Unfenced treatments are indicated by green and fenced treatments by red symbols and ellipse. Large symbols indicate the centroid of the plots. The ellipse indicates the 95% confidence interval of the treatment plots. The length of the blue solid arrows is proportional to its importance, and the angle between two arrows reflects the magnitude of the correlation between variables. Response variables were coded as: Soil_comp = soil compaction, Soil_pH = soil pH, Litter_depth = litter depth, Litter_quality = litter quality, Saping_div = sapling Shannon diversity, Sapling_dens = sapling density, Rodent_acti = rodent activity, Inv_div = invertebrate Shannon diversity, Inv_biomass = invertebrate biomass, Tea_k = tea decomposition, Pine_k = pine decomposition, Oak_k = oak decomposition, Birch_k = birch decomposition, N_min = nitrogen mineralization and PO_min = phosphate mineralization. Fix factors in dashed blue arrows are superposed and were coded as: Age = treatment age and Ung_abundance = ungulate abundance.

Full size image

In terms of browsing, broadleaf species (for example, Sorbus aucuparia, Rhamnus frangula, Quercus robur, Amelanchier lamarckii and Betula pendula) were highly browsed compared to conifer species (for example, Pinus sylvestris and Larix kaempferi, Figure 4). Concerning invertebrate composition, Hemiptera (aphids and plant-hoppers), Megadrilacea (earthworms), Diplopoda (millipedes), Hymenoptera (bees and ants), Coleoptera (beetles) and Geophilomorpha (centipedes) were associated with fenced plots and treatment age, whereas Diptera (flies) and Araneida (spiders) were associated with unfenced plots and with high ungulate abundance (Online Appendix A.2).

Figure 4
figure 4

Browsed stems result in percentage for all species across all unfenced plots in the Veluwe, the Netherlands. All stems above 15 cm in height and below 220 cm were used for this analysis. The number of replicate stems is showed in parenthesis. Browsing was scored as browsed or not browsed. Browsing intensity differed significantly among species (Kruskal–Wallis test, x2 = 432, d.f. = 12, p < 0.001). Some species (Quercus petrea, Betula pubescens, Pseudotsuga menziesii, Prunus serotine and Salix caprea) could not be quantified precisely because they had a low stem number.

Full size image

Fencing significantly affected nine out of 15 response variables, as indicated by the GLMM analyses (Figure 5, Table 1, Online Appendix A.3). Plots with ungulates (that is, unfenced) compared to plots without had significantly higher soil compaction (median = 375 vs. 294 kPa, Figure 5A), and lower litter depth (M = 2 vs. 4 cm, Figure 5B), sapling diversity (M = 0.71 vs. 1.05 H’, Figure 5C), sapling density (M = 10 vs. 35 ind m−2, Figure 5D), rodent activity (M = 1.5 vs. 4 bags removed, Figure 5E), invertebrate biomass (M = 1.53 vs. 1.65 log10mm3 kg−1, Figure 5F) and decomposition rate of tea bags (M = 0.01 vs. 0.02 k, Figure 5G), pine litter (M = 0.61 vs. 0.71 k, Figure 5H) and birch litter (M = 0.63 vs. 0.81 k, Figure 5I). All the other variables, such as phosphate and nitrogen mineralization rates, did not differ significantly between unfenced and fenced plots.

Figure 5
figure 5

Differences between treatment effect (fencing = 1 vs. unfencing = 0) for the nine significant response variables. All the individual data points (dots) are plotted as well. Significance (p < 0.05) was tested with general linear mixed models (GLMM, see Table 1). Rodent activity refers to the number of bags altered by rodents in the forest. Tea bag decomposition represents the combining values of green (fast) and red tea (slow) for one moment in time.

Full size image

Cascading Effects

There are several direct pathways through which wild ungulates triggered cascading effects in temperate forests (Figure 6, Online Appendix A.4). The PAM showed that in unfenced plots, ungulates significantly increased soil compaction (standardized regression coefficient β = 0.48) and decreased litter depth (β = − 0.44), sapling density (β = − 0.42) and diversity (β = − 0.51). Litter quality negatively predicted litter depth (β = − 0.27). These soil and vegetation variables had, in turn, relationships with other trophic levels. Sapling density positively predicted rodent activity (β = 0.42). Soil pH negatively predicted invertebrate biomass (β = − 0.45). Litter depth positively predicted invertebrate diversity (β = 0.33), whereas soil compaction (β = − 0.38) and soil pH (β = − 0.37) negatively predicted invertebrate diversity. Ungulates did not predict litter quality; similarly, tea bag, pine, oak and birch litter decomposition rate were not related to either invertebrate diversity or invertebrate biomass. Regarding nutrient mineralization, invertebrate biomass positively predicted nitrogen mineralization (β = 0.54).

Figure 6
figure 6

Cascading effects of wild ungulate presence on lower trophic levels, assessed with a fencing experiment of different ages. The results of a pathway analyses mode (PAM) are shown. The figure is based on the results of Online Appendix A.3. All model components were measured with N = 24. Solid arrows indicate significant (p < 0.05) relationships and dashed lines non-significant relationships. Numbers next to arrows are standardized regression coefficients, which are only given for significant relationships and inside the boxes the coefficient of determination (R2) are specified for the significant models only.

Full size image

Treatment Age and Ungulate Abundance

The PCA analysis (Figure 3) indicates that the superimposed variables: treatment age and ungulate abundance, had important associations with forest respond variables. For example, treatment age was positively associated with invertebrate biomass, nitrogen mineralization and phosphate mineralization and negatively associated with litter quality, oak decomposition and soil pH. Whereas, ungulate abundance was positively associated with soil compaction and negatively associated with sapling density, sapling diversity, rodent activity, litter depth and tea decomposition. The GLMM results indicate that treatment age did not have relationships with any of the response variables that directly interact with ungulates (Table 2), except for a slight non-significant relationship with litter depth (β = 1.60), although ungulate abundance had negative relationships with sapling diversity (β = − 0.01), sapling density (β = − 0.46) and litter depth (β = − 0.02).

Table 2 Results for the General Linear Mixed Models (GLMM) Between Ungulate Abundance and Treatment Age with Six of the Response Variables that Ungulates Have a Direct Effect
Full size table

Discussion

Ungulates are important ecosystem engineers that interact with forests at different levels. Yet the direction of the relationships is context dependent, meaning that ungulate abundance, primary productivity, food web complexity, apex predator presence and small consumer abundance can mediate the interactions between ungulates and forests (Forbes and others 2019). For our study area, we demonstrated that ungulate exclusion had multiple effects on forest ecosystems that appeared to cascade to different trophic levels. Ungulate exclusion changed the sapling composition, soil conditions and the invertebrate decomposer community, but did not have a significant effect on tree litter decomposition and nutrient mineralization rates. Further, in comparison with treatment age, ungulate abundance played a major role in the strength of the cascading effects, suggesting that ungulates are important ecosystem engineers, although in our system they had a limited role in nutrient cycling.

Fencing Effect

We found some support for our prediction that ungulate exclusion increases diversity on lower trophic levels and decomposition rate of litter, but also decreased soil compaction. Litter depth, sapling diversity, sapling density, rodent activity, invertebrate biomass, tea bag, pine and birch decomposition (Table 1, Figures 3, 5) were indeed significantly higher in the fenced compared to the unfenced plots. This implies that browsing and trampling directly and/or indirectly altered the vegetation and the composition of soil invertebrates (Gill and Beardall 2001; Allombert and others 2005b; Kuijper and others 2010). Invertebrate decomposers were mainly associated with fenced plots, whereas invertebrate predators were associated with unfenced plots (Online Appendix A.2). Soil compaction was significantly higher in unfenced plots, indicating that trampling exerted pressure on soil which directly translated to a higher soil density (Ramirez and others 2019). Fencing had a non-significant effect on N and PO4 soil mineralization rates, probably because nutrients tend to mineralize at initial stages of decomposition, and a three-month harvest might therefore be too long to still detect mineralization effects (Hasanuzzaman and Hossain 2014). Overall, the fencing experiment suggested that ungulates regulated forest succession into a poor state, characterized by a degraded trophic chain, high browsing incidence on palatable species and limited litter decomposition rates (Figures 3, 4, 5, Online Appendix A.2).

Above-ground Cascading Effects

We used PAM to analyze the pathways through which ungulates have cascading effects. Ungulates reduced understory stem density (Figure 6, Online Appendix A.4, β = − 0.42) and tree species diversity (β = −0.51) (Allombert and others 2005a, b). Regeneration density had, in turn, a positive relationship with rodent activity (β = 0.42), possibly because it protects rodents against extreme climatic conditions and provides refuge against predators (Flowerdew and Ellwood 2001). These results are confirmed by a study conducted in the UK, where bank vole (Myodes glareolus) and wood mice (Apodemus sylvaticus) population was lower in animal unfenced plots compared to fenced plots, and population size was directly related to the number of shrubs (Buesching and others 2011). Our results suggest that ungulates indirectly decreased rodent activity by shifting forest structure toward a less dense understory.

Below-ground Cascading Effects

We hypothesized that ungulates would affect soil macro-invertebrates by shifting tree species composition in the forest understory, which can reduce substrate quantity (litter depth) and soil quality (compaction). In sequence, our results showed that ungulate presence reduced invertebrate diversity by first reducing litter depth (Figure 6, β = − 0.44) and by increasing soil compaction (β = 0.48). In turn, litter depth (β = 0.33) increased invertebrate diversity, whereas soil compaction (β = − 0.38) decreased it. These results suggest that through trampling, ungulates directly reduced the litter layer depth by mixing litter with soil (Bruinderink and Hazebroek 1996; Hobbs 1996). The litter layer plays a vital role for soil invertebrates because it represents an important food source for decomposers and it also acts as a barrier that controls for humidity, temperature and light on the forest floor (Mills and Macdonald 2004). The negative relationship between soil compaction and invertebrate diversity agrees with the findings of Lavelle and others (2006), where invertebrate abundance decreased with increasing soil compaction and its side properties: limitations on soil water storage, soil aeration and invertebrate movement (Lal 1988; Althoff and Thien 2005).

The prediction that ungulates would reduce litter decomposition and mineralization by limiting the quantity and diversity of soil invertebrates was only partly supported (Figure 6, Online Appendix A.4). Invertebrate diversity did not have a significant relationship with decomposition and mineralization rates, whereas invertebrate biomass had only a positive relationship with nitrogen mineralization rate (β = 0.54). In agreement, an experimental study showed that invertebrate richness increased decomposition and mineralization rates (Jonsson and Malmqvist 2000), probably because of the cooperation among the decomposer guild (Gessner and others 2010). However, not all studies support the hypothesis of fast decomposition rate due to species cooperation (Hättenschwiler and others 2005), probably because environmental variables such as temperature and precipitation play a major role in litter decomposition and mineralization rates (Zhang and others 2008). In the case of this study, only surveying for macro-invertebrates could also had reduced the power of our results because decomposition is also driven by meso-fauna, fungi and bacteria (Aerts 2006).

Ungulate Abundance and Treatment Age

Since cascading effects are not only mediated by ungulate presence, we also tested for the relationship between ungulate abundance and treatment age with lower trophic levels. Litter depth was only weakly related to treatment age (Table 2, β = 1.60), suggesting that age might trigger cascading effects in this temperate forest. On the other hand, ungulate abundance significantly predicted sapling diversity (β = − 0.01), sapling density (β = − 0.46) and litter depth (β = − 0.02). The negative relationships between ungulate abundance and forest properties are widely discussed in the literature, and a meta-analysis recommends an ungulate density below 12 roe deer km−2 in order to conserve temperate forest structure, composition and functioning (Ramirez and others 2018). Considering that the average density in this forest site is 13 ungulates km−2, it is arguably to be expected that most of the interactions between ungulates and forest are negative.

Feedback Loops Promoted by Browsing

Feedback loops represent regulating forces on natural systems that make exploitative and symbiotic interactions persistent in time (Patten and Odum 1981), and the combination of different feedback loops allows systems to maintain their stability (Soto-Ortiz 2015). Feedback loops happen when the output of a particular ecosystem mechanism subdues the input of the same system. For example, ungulate browsing can trigger a positive feedback on savannah ecosystems, which leads to more intense fires by increasing the fuel load of grass and thus facilitating for grass biome (Van Langevelde and others 2003). Similarly, our results suggest that ungulates, through browsing on palatable species, may arrest forest ecosystems into an early successional stage with low primary productivity, as in Ramirez and others (2019). For instance, by browsing on palatable species (Figure 4), the quantity and quality of litter production may be reduced leading to a decrease in soil invertebrate diversity and composition (Figure 3, Online Appendix A.2), which in turn could theoretically stall decomposition and mineralization rates. A slow nutrient cycle in the soil, results, in turn, in less nutrients for plant development and selects for plant species with low nutrient requirements and nutrient-poor litter (Figure 1) (Augustine and McNaughton 1998).

A fair amount of literature has addressed ungulate effects on multi-trophic systems across different biomes: from boreal forests to savannahs. The majority of studies agree that ungulates exert through browsing and trampling a strong top-down control on vegetation, which triggers negative cascading effects on lower trophic levels and ecosystem processes. These effects promote slow feedback loops that reinforce the establishment of less palatable species (for example, conifer) that require a limited amount of soil nutrients (Pastor and others 1993; Ritchie and others 1998). Yet, a group of studies demonstrated that under certain environmental conditions, ungulates could trigger positive cascading effects on systems, which can lead to fast feedback loops (Pastor and others 1988). For example, Frank and others (2000) demonstrated that ungulates on unfenced compared to fenced plots in Yellowstone National Park, promoted nitrogen availability by stimulating microbial activity and microbial turnover rates. Similarly, through a literature review Augustine and McNaughton (1998) proposed with empirical evidence that ungulate browsing may also establish the dominance of highly palatable species, especially when there is a high level of nutrient inputs which facilitate individuals to recover from tissue loss or by when there is a good balance between tolerant and intolerant tree species to browsing. Finally, Pringle and others (2007) revealed that ungulate presence decreased the abundance of trees, lizards and arthropods in a series of ungulate fenced compared to unfenced plots that ranged in productivity in an African savannah; yet, the strength of these trophic effects decreased with increasing primary productivity.

The growing amount of evidence across biomes suggests that external factors such as climate and primary productivity may be more important on mediating trophic interactions and ecosystem processes (Hobbs 1996). Considering that our project was conducted in a nutrient deficient system and with limited primary productivity, we expected and found that ungulate presence had strong cascading effects on lower trophic levels, yet we only found weak relationships with nutrient mineralization and none with litter decomposition. Suggesting that not only external factors such as climate and primary productivity mediate trophic interactions, but actually, the community organization and structure of each trophic level may also have an overriding force.

Conclusions

Our findings suggest that wild ungulates have strong direct and indirect cascading effects on forest ecosystems. By browsing on palatable tree species, ungulates homogenize lower trophic levels, which have the potential to trigger slow feedback loops in this system. This can eventually arrest forest succession to an early stage mostly composed by light demanding tree species—which are less palatable. Surprisingly, nutrient cycling was not affected by ungulates, implying that forest functions are highly resilient to browsing and trampling. These cascading effects are amplified with increasing ungulate abundance, and thus, limiting ungulate populations might improve temperate forest resilience.