Ungulate Species and Abundance as well as Environmental Factors Determine the Probability of Terminal Shoot Browsing on Temperate Forest Trees

: Ungulate browsing is a major factor inﬂuencing tree regeneration. However, it is unclear if the observed increase in ungulate abundance in Central Europe implies increased browsing, and which other factors inﬂuence the incidence of browsing. We investigated the impact of forty variables (site, climate, forest and ungulates) on the probability of leader shoot browsing of six tree species which are frequent in Switzerland. The analysis was based on a large dataset including 49 monitoring areas, each containing 25–64 circular plots, in which 10 to 130 cm tall seedlings were repeatedly assessed. Browsing probability was estimated for each plot and year by mixed e ﬀ ects logistic regression and used as a response in random forests to disentangle the inﬂuence of the explanatory variables. Browsing probability was positively correlated with ungulate density measures (number culled by hunting or found dead) for all six tree species. Where beyond roe deer, some red deer and / or chamois were present, the browsing probability was higher. Small timber tree stands had less browsing than young growth and thicket stands. Seedlings tended to be more frequently browsed in stands with > 80% canopy shading. Browsing increased with increasing understory cover, independent of vegetation category. In conclusion, browsing is a multifactorial phenomenon and ungulate density estimates alone do not explain the whole browsing probability.


Introduction
Ungulates such as red deer (Cervus elaphus L.), roe deer (Capreolus capreolus L.) and chamois (Rupicapra rupicapra L.) strongly affect forest regeneration in many temperate and boreal forests (e.g., References [1,2]) since woody plants often make up a large part of their diet [3][4][5]. In central Europe, ungulate abundance has increased in recent decades and is likely to continue to do so in the future [6].
Deer and chamois are selective browsers, or at least intermediate feeders [7]. Thus, they select tree seedlings by actively choosing certain tree species [8]. The preferred ones are used disproportionately to their availability [8]. Among the coniferous tree genera, Abies are typically the most preferred in northern temperate forests of Europe, while Picea and Pinus are usually least preferred [9]. In Switzerland, annual leader shoot browsing, i.e., feeding on the apical bud or terminal shoots of tree seedlings,

Study Sites and Sampling Design
In order to monitor the frequency of ungulate browsing on tree regeneration, leader shoot browsing has been repeatedly (i.e., annually or every second year) recorded on more than 200 indicator areas (IA) throughout Switzerland [38]. An IA typically comprises a continuous forest area with different developmental stages and forest types, which is characteristic for the respective ungulate management unit [39]. In this study, we exploit this large dataset by selecting 49 IA which (i) have been monitored for at least 10 years, (ii) have been recently monitored (2008 or later), (iii) do not have any known methodological deficits such as translocation of sampling plots and (iv) contained at least 10 sampling plots (at any point in time) with presence of Abies seedlings, among other species, in order to be able to compare browsing of the different species. The selected IA were all located in north-eastern and central Switzerland in 5 cantons (Figure 1), i.e., 36 IA in St. Gallen (SG), 10 in Thurgau (TG) and one each in Nidwalden, Obwalden and Schwyz.
Forests 2020, 11, x FOR PEER REVIEW 3 of 20 In order to monitor the frequency of ungulate browsing on tree regeneration, leader shoot browsing has been repeatedly (i.e., annually or every second year) recorded on more than 200 indicator areas (IA) throughout Switzerland [38]. An IA typically comprises a continuous forest area with different developmental stages and forest types, which is characteristic for the respective ungulate management unit [39]. In this study, we exploit this large dataset by selecting 49 IA which (i) have been monitored for at least 10 years, (ii) have been recently monitored (2008 or later), (iii) do not have any known methodological deficits such as translocation of sampling plots and (iv) contained at least 10 sampling plots (at any point in time) with presence of Abies seedlings, among other species, in order to be able to compare browsing of the different species. The selected IA were all located in north-eastern and central Switzerland in 5 cantons (Figure 1), i.e., 36 IA in St. Gallen (SG), 10 in Thurgau (TG) and one each in Nidwalden, Obwalden and Schwyz. In each IA, 25-64 circular sampling plots were established on a rectangular grid with a fixed distance of 100 to 200 m (Figure 1), and the plot centers were permanently marked. This resulted in a total of 1704 plots. The plot radius was either 2.0 or 5.0 m, depending on seedling abundance in an IA, and corrected for slope ( Figure 1, [40]).

Long-Term Field Data Collection
Each plot was sampled yearly or every other year in late winter or early spring before budburst for 10 to 16 years, yielding 6 to 10 inventories in each plot. Tree seedlings between 10 and 130 cm were counted in four height classes (hc) (hc1: 10-40 cm, hc2: 41-70 cm, hc3: 71-100 cm and hc4: 101- In each IA, 25-64 circular sampling plots were established on a rectangular grid with a fixed distance of 100 to 200 m (Figure 1), and the plot centers were permanently marked. This resulted in a total of 1704 plots. The plot radius was either 2.0 or 5.0 m, depending on seedling abundance in an IA, and corrected for slope ( Figure 1, [40]).

Long-Term Field Data Collection
Each plot was sampled yearly or every other year in late winter or early spring before budburst for 10 to 16 years, yielding 6 to 10 inventories in each plot. Tree seedlings between 10 and 130 cm were counted in four height classes (hc) (hc1: 10-40 cm, hc2: 41-70 cm, hc3: 71-100 cm and hc4: 101-130 cm), two browsing classes (leader shoot browsed or not browsed in the last year) and nine species categories. The species categories were Abies (A. alba Mill.), Picea (P. abies (L.) H.Karst.), other coniferous species, Fraxinus (F. excelsior L.), Acer (containing A. pseudoplatanus L., A. platanoides L. and rarely A. campestre L.), Fagus (F. sylvatica L.), Quercus (containing Q. petraea (Mattuschka) Liebl. and Q. robur L.), Sorbus (containing S. aucuparia and rarely S. aria (L.) Crantz) and other broad-leaved tree species. Quercus was only assessed in IA in the canton TG and Sorbus in all other cantons. As the categories 'other broad-leaved trees' and 'other coniferous trees' contained different species, they were omitted in species-specific analyses. The assessment of seedlings always started in the north direction and proceeded clockwise around the center of the sampling plot. If the number of seedlings in a plot exceeded 30, counting was stopped to save time and the azimuth was recorded ( Figure 1). Thus, an observation in a certain year and at a specific location (IA:plot) consisted of the number of browsed and unbrowsed seedlings per height class and species category, as well as a radius and azimuth value.

Additional Assessment of Variables
In order to assess factors that may influence browsing, we complemented the long-term monitoring data with variables related to location, climate, forest characteristics and ungulates (Appendix A Table A1). All these variables, except ungulate-related variables, were determined at only one point in time (e.g., forest type) or as an average over time (e.g., precipitation) and were assumed to be constant over time for a specific sampling location (i.e., at IA:plot).
Location variables were derived from a geographic information system (GIS) using swisstopo data DHM25 (with a 25 m grid, [41]). Aspect was converted to a factor with four levels according to the cardinal directions. Daily climate data of the period 1990 to 2010 provided by Meteoswiss from all available climate stations was interpolated over the 100 m Digital Elevation Model of Switzerland with Daymet according to Thornton et al. [42] to get means for annual precipitation, precipitation from April to June, annual temperature and degree day sum for each plot. Degree day sums were calculated [43] with a minimum threshold of 5.56 • C.
Forest characteristics for each plot were retrieved from geographic information systems of the cantonal authorities. Maps of forest associations were used to estimate soil humidity and acidity according to the ecograms in Frehner et al. [44]. Abies-Fagus associations were classified as coniferous forests, due to their aspect in winter.
Ungulate densities in each IA were approximated with the number of roe and red deer and chamois culled by hunting or found dead for other reasons (such as traffic accidents, wounded after being shot, agricultural machinery, predators, or natural deaths like crashes, diseases, old age) per year and forest area within each hunting district (data provided by cantonal authorities). As it is not possible to identify the ungulate species responsible for leader shoot browsing in such a large repeated assessment, we looked for an ungulate density index (UDI) to represent the overall impact of all ungulate species present. This index should take into account the body mass but also the ratio of woody vegetation to total diet [31]. Following Motta [31], we calculated UDI = 1 5 roe deer + 1 4 chamois + red deer, where all densities refer to 100 ha. In addition, the ungulate composition variable quantified whether and how many red deer and chamois were present in comparison to roe deer (0%, <5%, ≥5%; Appendix A Table A1).

Supplementary Field Data Collection in 2012
To answer the research questions (iv)-(vi), a subsample of 660 plots in 40 IA was revisited in 2012 to record light and vegetation cover (Appendix A Table A1). Due to the importance of Abies, plots were only selected if they contained at least three Abies seedlings at any point in time. Furthermore, the plots had to be in the cantons SG and TG. Some plots were irretrievable due to missing terrain marks or difficult terrain conditions after a local windthrow event.
Shading at each plot was recorded at 40 and 130 cm above ground as the proportion of the sky covered by both the canopy after complete leaf opening and the local horizon, corresponding to the proportion of black in a black and white fisheye picture taken vertically above the plot center (i.e., criterion "Shading" in the Swiss forest inventory, cf. Reference [45]).
The canopy above a plot was optically separated into eight sectors and each of these sectors was recorded as gap if less than half was covered by treetops (Supplementary Figure S1). The variable 'gap size' is the sum of the gap sectors for each plot, i.e., zero if no opening is present in the canopy, and higher in larger gaps. Different weighting of the gap sectors was used to approximate snow depth and light transmittance. The latter was based on calculations of total transmitted light (i.e., the sum of direct and indirect transmitted light) for openings of single sectors with different slope and aspect with the Gap Light Analyzer software [46].
Basal areas of living trees, dead trees (snags and fresh lying logs) and stumps were quantified using angle count sampling with counting factor 4 [47].
Total vegetation cover, as well as the cover of blackberry, raspberry, ferns, graminoids, perennial herbs, blueberry and shrubs, were estimated within a plot radius of 2 m. All these vegetation covers and seedling cover were recorded on a scale introduced by Londo [48]. Based on this data, no statement can be made on species richness in the plots but it can tell whether ground vegetation cover affects browsing on tree seedlings.

Data Analysis
We used a two-step process to analyze the extensive data. Browsing probabilities (PB) were first estimated for each species, plot and year using a mixed effects logistic regression with proximal explanatory variables (such as species, year and height class), and accounting for the binary nature of the outcome (seedling browsed/not browsed). These PB estimates were then used as the response variables in random forest models with a large number of more distal explanatory variables (Appendix A Table A1). The two-step process combined the benefits of each method by using the mixed effects logistic regression to model the temporal and spatial correlation within the data and random forests to assess the predictive effects of a large number of explanatory variables.
Mixed effects logistic regression allowed for estimating PB for each plot, including those containing few seedlings where the inventory yielded highly stochastic browsing percentages (e.g., on plots with one seedling only, which is either browsed or not). We included fixed effects for species, height class and year (without intercept), and random intercepts for IA, IA:plot, species:IA:plot, year:IA:plot, hc:IA:plot, species:year, species:hc and year:hc. The random effect structure was chosen to model potential correlations, allowing for, e.g., species-or year-specific variations of the browsing level at a given location. It introduces rudimentary spatial and temporal correlations via the location-and year-specific random intercepts. The three-way interactions were introduced to model correlations within a specific year, species and height class at a given location (defined by plot within an IA). They were based on the assumptions that (i) seedlings at a given location of a specific species may share a common browsing level (e.g., some species could be more selected due to microsite conditions, ground vegetation or species composition), (ii) seedlings at a given location in a given year may share a common browsing level (e.g., inter-annual climatic variability could lead to small-scale shifts in the spatial use of habitats), and/or (iii) seedlings at a given location with a given height may share a common browsing level (e.g., a recurrent snow cover could affect PB of a certain height class but not the species browsed). The two-way interactions, species:year, species:hc and year:hc, were introduced to capture further correlations within species in a year or height class or between years and height class. This allowed for including the possibility that, e.g., one species (but not the others) might have been generally heavily browsed in one year at all locations. Relationships between IAs were not included, i.e., the 49 IAs were assumed to be independent.
We did not directly further add more distal explanatory variables to the mixed effects logistic regression to avoid a complicated and potentially biased variable selection process and the restriction to linear effects. Instead, we used random forests, a specialized regression method based on an ensemble of regression trees. Random forests do not make any assumptions about the distribution of the response or the effect of the explanatory variables on the response and can cope with many correlated explanatory variables [49]. Separate random forests were fitted for Picea, Abies, Fagus, Fraxinus, Acer and Sorbus (but not for Quercus, which were only present in the canton TG). We used year, height class and all variables shown in Appendix A Table A1 as explanatory variables, restricting the analysis to the 40 IAs and 660 plots where supplementary field data had been collected in 2012. To ensure that this restriction did not have a major effect on the results, additional random forests models were fitted in a sensitivity analysis, using only the subset of explanatory variables available for all 49 IAs (i.e., excluding the variables of the supplementary field data collection in 2012, cf. lower part of Appendix A Table A1).
To quantify the variation within the data, random forests were fitted to resampled datasets comprising only 25 of the 40 IAs. The predictive importance of a variable on the response was quantified using the variable importance calculated according to Breiman's method (also called permutation importance, [50]). Marginal effects were assessed via the partial dependence function (function partialPlot of the R package randomForest, see below) and presented in partial dependence plots. This procedure was extended to more than one variable, leading to their joint marginal effects. Separate subset analyses were conducted for the three ungulate composition categories.
A sensitivity analysis was conducted to account for the uncertainty in estimating PB. The probability of an observation to be included in an individual regression tree of the random forest models was weighted by the inverse of its standard error from the mixed effects logistic regression. That is, more precise observations were selected more often and got more weight in the subsequent random forests.
We performed all statistical analyses with R version 3.0.2 [51]. The logistic mixed model was fitted using the function glmer with the family specification 'binomial' (link = 'logit') of the R package lme4, version 1.0-5 [52]. Random forests were fitted with the R package randomForest, version 4.6-7 [53], using 500 'trees' and a minimal node size of 5.

Distribution of Ungulates and Browsing Probability
The annual number of dead animals, as an indicator of abundance, varied between 0 and 27.5 km −2 in roe deer, 8 km −2 in chamois and 4 km −2 in red deer (Appendix A Table A1 and Supplementary Figure S2). Red deer and chamois were both absent in 35% of the IA. Roe deer was the only ungulate species present in 9 out of 10 IA in the canton TG and in 8 out of 36 IA in the canton SG. These IA were predominantly located in the Swiss Plateau.
In total, 275,866 seedlings (between 11,577 and 192,252 per height class) on 1704 plots were assessed for leader shoot browsing (Appendix A Table A2). Depending on tree species and height class, between 4% and 67% of the tree seedlings were browsed. Both modelled PB and observed percentage of browsing are estimates for the unknown browsing probability. If many seedlings were present, both estimates were very similar, but if only few were present, they differed considerably (Figure 2, e.g., situations with an observed proportion of 0, 0.5 or 1, which often occurred if zero, one or two seedlings were present). In such cases, PB took into account information from other years, height classes and also from plots in the same IA and other species (see estimated model parameters in Appendix A Table A3), whereas the observed percentage of browsing did not. Figure S2). Red deer and chamois were both absent in 35% of the IA. Roe deer was the only ungulate species present in 9 out of 10 IA in the canton TG and in 8 out of 36 IA in the canton SG. These IA were predominantly located in the Swiss Plateau.
In total, 275,866 seedlings (between 11,577 and 192,252 per height class) on 1704 plots were assessed for leader shoot browsing (Appendix Table A2). Depending on tree species and height class, between 4% and 67% of the tree seedlings were browsed. Both modelled PB and observed percentage of browsing are estimates for the unknown browsing probability. If many seedlings were present, both estimates were very similar, but if only few were present, they differed considerably (Figure 2, e.g., situations with an observed proportion of 0, 0.5 or 1, which often occurred if zero, one or two seedlings were present). In such cases, PB took into account information from other years, height classes and also from plots in the same IA and other species (see estimated model parameters in Appendix Table A3), whereas the observed percentage of browsing did not. The estimated PB for Fagus and Picea tended to be slightly lower than the observed browsing probabilities ( Figure 2). However, most of the plots with seedlings of these two species are close to The estimated PB for Fagus and Picea tended to be slightly lower than the observed browsing probabilities ( Figure 2). However, most of the plots with seedlings of these two species are close to the lower left corner of Figure 2 and the few strata with high observed browsing probabilities only weakly contribute to the model fit.
The observed and modelled PB highly varied between IA, and both IA with rare browsing (e.g., IA 12, 45 and 49) and frequent browsing (e.g., IA 18, 30 and 42) were observed ( Figure 3, left panels). Nevertheless, the species-specificity was remarkably constant: Sorbus and Quercus were most often browsed, followed by Acer and Fraxinus, then by Abies and finally Fagus and Picea (Figure 3). This was also reflected by the odds ratios from the mixed effects logistic regression (Appendix A Table A3). Seedlings of small (hc1, 10-40 cm) and especially of tall height classes (hc4, 101-130 cm) were less often browsed than medium-sized ones (hc2, 41-70 cm and hc3, 71-100 cm; Figure 3, right panels). Accordingly, the odds for browsing were higher in hc2 and hc3 than in hc1, but smaller in hc4 (Appendix A Table A3).

Variables Influencing Browsing Probability
For the six tree species, a similar set of variables was important for predicting PB (Figure 4), and their marginal effects on PB showed common patterns ( Figure 5 for selected variables for Abies, Supplementary  Figures S3-S9 for all variables and species). For the rarely browsed species Picea and Fagus (Figure 3), all variables showed a lower variable importance compared with Abies and the other more frequently browsed deciduous species (Figure 4). There were also smaller differences in the marginal effects on PB (Supplementary Figures S3-S9). This can be explained by the fewer strata with high browsing probabilities (cf. Figure 2).
The results were confirmed by both sensitivity analyses. The analysis conducted with the 18 explanatory variables available for all 49 IAs (cf. Appendix A Table A1, but with UDI excluded) showed similar patterns for predictive importance (Supplementary Figure S10) and marginal effects (not shown). The random forest models that accounted for the uncertainty of PB-by including weights relative to the precision of the PB estimates of the mixed effects logistic regression-yielded almost the same variable importance (Supplementary Figure S11) and marginal effects (Supplementary Figures S12 and S13) as the main model. Nevertheless, the species-specificity was remarkably constant: Sorbus and Quercus were most often browsed, followed by Acer and Fraxinus, then by Abies and finally Fagus and Picea (Figure 3). This was also reflected by the odds ratios from the mixed effects logistic regression (Appendix Table A3). Seedlings of small (hc1, 10-40 cm) and especially of tall height classes (hc4, 101-130 cm) were less often browsed than medium-sized ones (hc2, 41-70 cm and hc3, 71-100 cm; Figure 3, right panels). Accordingly, the odds for browsing were higher in hc2 and hc3 than in hc1, but smaller in hc4 (Appendix Table A3).

Variables Influencing Browsing Probability
For the six tree species, a similar set of variables was important for predicting PB (Figure 4), and their marginal effects on PB showed common patterns ( Figure 5 for selected variables for Abies, Supplementary Figures S3-S9 for all variables and species). For the rarely browsed species Picea and Fagus (Figure 3), all variables showed a lower variable importance compared with Abies and the other more frequently browsed deciduous species (Figure 4). There were also smaller differences in the marginal effects on PB (Supplementary Figures S3-S9). This can be explained by the fewer strata with high browsing probabilities (cf. Figure 2).   The results were confirmed by both sensitivity analyses. The analysis conducted with the 18 explanatory variables available for all 49 IAs (cf. Appendix Table A1, but with UDI excluded) showed similar patterns for predictive importance (Supplementary Figure S10) and marginal effects (not shown). The random forest models that accounted for the uncertainty of PB-by including weights relative to the precision of the PB estimates of the mixed effects logistic regression-yielded almost the same variable importance (Supplementary Figure S11) and marginal effects ( Supplementary  Figures S12 and S13) as the main model.

Is Browsing on Tree Seedlings Positively Correlated with the Density of Ungulates?
Among the distal variables, the ungulate density index (UDI) and the approximated roe deer density were important predictors ( Figure 4) and had a positive effect on PB for all species ( Figure 5 and Supplementary Figure S6). PB increased almost linearly within the common range of approximated roe deer densities (between the first and ninth decile, Figure 5 and Supplementary Figure S6). The effect of chamois and of red deer was similar but weaker. The effect of UDI was similar for all height classes, as indicated by the joint marginal effects of UDI and height class (Supplementary Figure S14). Roe deer was also among the most important variables in all random forests conducted with the 18 explanatory variables (cf. Appendix Table A1, but with UDI excluded) available for all 49 IAs (Supplementary Figure S10).

Does Browsing on Tree Seedlings Depend on Ungulate Species?
In contrast to UDI, ungulate composition was not among the most important variables ( Figure  4 and Supplementary Figure S10 for all 49 IAs), but PB was always smallest if only roe deer was present (category 0 in Figure 5 and Supplementary Figure S6). Subset analysis for observations with roe deer only (category 0) compared to observations with some additional red deer and chamois (categories 1 and 2 for less or more than 5% of roe deer, respectively) revealed almost the same marginal effects, just at lower PB values for roe deer only (Supplementary Figure S15). The difference between only roe deer and few other animals reached ≥10% of PB for Abies, Acer, Fraxinus and Sorbus (Supplementary Figure S15). Furthermore, UDI had slightly higher power to explain variation in PB than roe deer density for almost all tree species (Figure 4) and PB increased more strongly in dependence of UDI than of roe deer density ( Figure 5 and Supplementary Figure S6). The proximal variables height class and year had the largest or second largest predictive importance for five of the six species. Medium height classes were clearly associated with higher PB ( Figure 5 and Supplementary Figure S3). The effect of year may be seen as a surrogate for unobserved time-dependent variation.

Is Browsing on Tree Seedlings Positively Correlated with the Density of Ungulates?
Among the distal variables, the ungulate density index (UDI) and the approximated roe deer density were important predictors ( Figure 4) and had a positive effect on PB for all species ( Figure 5 and Supplementary Figure S6). PB increased almost linearly within the common range of approximated roe deer densities (between the first and ninth decile, Figure 5 and Supplementary Figure S6). The effect of chamois and of red deer was similar but weaker. The effect of UDI was similar for all height classes, as indicated by the joint marginal effects of UDI and height class (Supplementary Figure S14). Roe deer was also among the most important variables in all random forests conducted with the 18 explanatory variables (cf. Appendix A Table A1, but with UDI excluded) available for all 49 IAs (Supplementary Figure S10).

Does Browsing on Tree Seedlings Depend on Ungulate Species?
In contrast to UDI, ungulate composition was not among the most important variables (Figure 4 and Supplementary Figure S10 for all 49 IAs), but PB was always smallest if only roe deer was present (category 0 in Figure 5 and Supplementary Figure S6). Subset analysis for observations with roe deer only (category 0) compared to observations with some additional red deer and chamois (categories 1 and 2 for less or more than 5% of roe deer, respectively) revealed almost the same marginal effects, just at lower PB values for roe deer only (Supplementary Figure S15). The difference between only roe deer and few other animals reached ≥10% of PB for Abies, Acer, Fraxinus and Sorbus (Supplementary Figure S15). Furthermore, UDI had slightly higher power to explain variation in PB than roe deer density for almost all tree species ( Figure 4) and PB increased more strongly in dependence of UDI than of roe deer density ( Figure 5 and Supplementary Figure S6).

Is Browsing on Tree Seedlings Increased in Early Developmental Stages?
Developmental stage was an important predictor ( Figure 4 and Supplementary Figure S10 for all 49 IA). PB tended to be highest in the young growth and thicket stage and lowest in small timber tree or multi-layered stands for all species (cf. developmental stage based on aerial imagery, Supplementary Figure S5, or assessed in the field data collection 2012, Figure 5 and Supplementary Figure S7). Stands dominated by medium and large timber trees, which had higher basal area (36.9 ± 12.1 m 2 , mean ± standard error) than stands in other developmental stages (25.2 ± 13.2 m 2 ), had higher PB than stands of small timber trees in most species. PB often reached a minimum at an intermediate basal area (Figure 5), except in Acer, where higher basal areas were related to lower PB (Supplementary Figure S7).

Does Browsing on Tree Seedlings Differ between Gaps and Closed Canopy Stands?
Light variables were not among the most important variables ( Figure 4) and no marked effects could be observed, in particular none related to gap size ( Figure 5 and Supplementary Figure S7). The joint marginal effects of gap size and seedling height class did not provide evidence for an interaction, indicating that this effect was similar for all height classes. Only for the shade-or moderately shade-tolerant species Abies, Picea, Acer and Fagus [54], did PB tend to increase in closed stands (>80% canopy shading, variable 'shading at 130 cm', Figure 5 and Supplementary Figure S7).

Does the Presence of Stumps, Snags and Lying Logs Enhance the Probability of a Seedling to be Browsed?
Basal area of stumps and dead trees (snags and lying logs) was not among the most important variables (Figure 4). PB tended to increase with both the number of stumps and dead trees for Abies and Picea, whereas no change or rather a small decrease was observed for the other species (Supplementary Figure S8). However, the threshold-like increase for Abies (and partly for Picea and Acer) at high basal areas of stumps and dead trees was only observed for a few plots (much less than 10% of the data and not in all 25 resamples) and may have been caused by windthrow or logging.

Do the Cover of Tree Seedlings and Ground Vegetation Affect Browsing on Tree Seedlings?
Important variables explaining variation in PB were seedling and vegetation cover (Figure 4), although the effect size was very small ( Figure 5, Supplementary Figures S8 and S9). Increasing seedling cover seemed to be related to decreasing PB at low densities and with the opposite or only small effects thereafter ( Supplementary Figures S5 and S8). Vegetation cover was associated with increasing PB in all tree species (Figure 5 and Supplementary Figure S9), a trend also observed for the cover by individual vegetation types (e.g., blackberry, raspberry, ferns, graminoids, perennial herbs, blueberry and shrubs, Supplementary Figure S9).

Browsing Is a Multifactorial Phenomenon
Using a large dataset from 49 forest areas (indicator areas) assessed over 10-16 years, we studied the simultaneous effects of many variables on the probability of ungulate browsing on seedlings of six tree species in Central Europe. Browsing turned out to be a multifactorial phenomenon. None of the explanatory variables alone had a large effect on browsing. The browsing probability varied by a few percent only if one specific variable changed.
However, two sensitivity analyses demonstrated the robustness of the results: the reduction in the dataset caused by the incomplete field data collection in 2012 did not introduce a detectable bias (Supplementary Figure S10) and neglecting the uncertainty in PB hardly influenced the results of the random forest models (Supplementary Figures S11-S13).
In our discussion, we will focus on the six research questions, although our models showed that further factors also influence PB, for example, the proximate driver's height class (with medium-sized seedlings being more frequently browsed) and year. As tree seedlings of all species will need to grow through all height classes to survive, we think that a discussion on at which stage they are affected by browsing is not of major interest. The large influence of the factor 'year' probably results from the fact that most of our variables were not available on a yearly basis (e.g., local climate variables), meaning that year is very likely a surrogate for yearly changes not captured by any covariate.

Is Browsing on Tree Seedlings Positively Correlated with the Density of Ungulates?
Our analysis shows browsing to be more frequent if ungulate densities (approximated with the number of roe, red deer and chamois culled by hunting or found dead for other reasons) are higher. Likewise, the incidence of leader browsing by deer was positively associated with an index of deer density at the site scale [55] or with deer presence at the plot scale [15,30]. In addition, Hothorn et al. [56] showed that the risk of deer-vehicle collisions is positively correlated to browsing intensity and to harvest numbers of ungulates.
The ungulate density variables and PB were within the frequently observed densities linearly related. Just at very high roe deer densities (more than 20 deaths km −2 ), a steep increase of browsing was observed (see Abies in Figure 5). This might reflect an increased use of tree seedlings due to competition for limited resources if roe deer is abundant. However, high densities were scarce (less than 10%) in our data and the uncertainty of the marginal effects was large.
It has to be kept in mind that browsing is species-specific. A certain number of ungulates (in our study, for example, ca. 18 dead roe deer per km 2 ) still cause occasional browsing on Picea or Fagus seedlings but can already lead to browsing of every fourth Abies and/or every third Sorbus seedling. This strong species preference may explain why in some studies with high ungulate levels, a hunting increase does not reduce browsing on the preferred species [36]. The relationships between ungulate abundance and browsing on different tree species are certainly not running in parallel.

Does Browsing on Tree Seedlings Depend on Ungulate Species?
We found higher browsing if not only roe deer but also red deer and/or chamois were present. This suggests that not only the density but also the composition of the ungulates influences the level of browsing. Onderscheka et al. [57] found more remains of coniferous trees in the stomach of red deer and chamois than in that of roe deer in the Principality of Liechtenstein, which borders our study area. In a study in the Swiss Prealps, red deer and chamois were filmed while browsing on leader shoots of tree seedlings in winter quite often, but roe deer only occasionally [58]. This was induced by seasonal shifts, as roe deer moved to lower elevations over the winter. These findings suggest that at least at certain sites, red deer and chamois browse more frequently on leader shoots than roe deer. However, ungulate diets vary largely between forest sites [4,22], and interactions between ungulate species that co-occur in an area can alter their browsing behavior [13]. We cannot exclude the alternative explanation that the roe deer in our indicator areas with multiple ungulate species browse more often on tree seedlings than do roe deer without competitors. That is, diverse ungulate assemblages could lead to more browsing on low-quality resources.
Browsing and ungulate density were more closely related for roe deer than for red deer or chamois. Indeed, roe deer occurred in all our indicator areas, and their densities covered a wider range. Thus, the roe deer effect might just have been easier to detect. More specialized studies based on DNA marker technology are needed to provide insight into which ungulate species is really browsing on tree seedlings [59,60]. This can facilitate effective targeting of management interventions when dealing with multiple ungulate communities [61].

Is Browsing on Tree Seedlings Increased in Early Developmental Stages?
Browsing was more frequent in young growth and thicket stands than in those dominated by small timber trees. However, the browsing probability was also high in stands dominated by medium or large timber trees. This may also explain why PB often reached a minimum at an intermediate basal area. We therefore have no clear evidence to support high browsing probabilities in early developmental stages. Young growth stands and thickets may not be particularly prone to browsing, whereas stands dominated by small timber trees often have little regeneration and are unattractive for ungulates.
Multi-layered stands were rare in our data but tended to have low browsing probabilities. Yet, multi-layered continuous forests are attractive habitats for ungulates [62] and are supposed to be more predisposed to silvicultural damage caused by browsing compared to other silvicultural systems [63]. Low regeneration density with a high proportion of Abies (which is more frequently browsed than Picea), rather than excessive browsing per se, might result in the supposedly high browsing damage in such stands (see also Reference [64]). Indeed, we found positive examples in Switzerland, where natural regeneration of Abies could grow up again well through ungulate browsing, almost exclusively in such multi-layered stands [65].

Does Browsing on Tree Seedlings Differ between Gaps and Closed Canopy Stands?
We did not find a clear and important association of browsing and canopy characteristics, in particular not for gap size, but seedlings in dark stands with canopy openness < 15-20% tended to be browsed more frequently. This is in agreement with studies where browsing was similar in logged windthrow gaps and adjacent forest stands [66] and along the whole gradient in between [21], or with a study in which more ungulate browsing occurred in very dark stands and more mice browsing in open conditions [34]. Similarly, but in a more pronounced way, ungulate browsing decreased with an increase in direct light [67] or was less frequent in uncleared windthrow gaps than in adjacent forest stands [26,68]. However, the relationship between PB and light is a matter of debate, as also the contrary was found. For example, the visitation frequency of ungulates was higher in gaps than in closed mature stands [69], and seedlings in gaps were therefore more frequently browsed [25]. Ungulates may select gaps if stands with high visual protection are located nearby [70], but avoid them in case of too much snow [34,71]. While our study is based on a large set of longitudinal data, detecting gap effects on browsing was probably difficult since they are often overruled by other factors such as surrounding stands or snow cover [72].

Does the Presence of Stumps, Snags and Lying Logs Enhance the Probability of a Seedling to Be Browsed?
The amount of fresh coarse woody debris (CWD; measured as basal area of snags, stumps and lying logs) neither had an important nor the same effect on PB of all tree species. Both coniferous species tended to be more frequently browsed with higher amounts of CWD. This trend is in line with a montane Picea abies snag stand, where earlier snow melt underneath or close to logs lead to seedlings in such microsites being browsed more often than seedlings distant to logs [29]. Similarly, in another study, the presence of woody debris increased deer browsing and reduced the abundance of Abies alba and Quercus spp. regeneration [28].
We found no clear interaction between snags, logs or stumps and the browsing probability in deciduous tree species. This is in agreement to findings in the Białowieża National Park outside the core area of wolves [73]. Inside the wolf core area, browsing intensity was reduced if the amount of CWD increased, suggesting that CWD can act as an 'ungulate escape impediment' [73]. Generally, CWD may act less as an obstacle for chamois and roe deer than for red deer [68]. Thus, the effect of logs and stumps on browsing seems to depend on several factors, including snow melt, predator presence and ungulate species, and definitively needs further research attention.

Do the Cover of Seedlings and Ground Vegetation Affect Browsing on Tree Seedlings?
Ungulate browsing tended to be more frequent if alternative food was available, i.e., if vegetation cover was higher. Browsing is generally assumed to be particularly high if palatable understory vegetation, like blackberry and raspberry, or seedlings of selected tree species are present [17,20]. However, we found the same pattern for all seven analyzed vegetation categories, although the importance ranks differed between tree species and vegetation categories.
Generally, seedlings are most easily detected and browsed in short ground vegetation, i.e., if they protrude over it [26,74]. Therefore, a dense and tall shrub layer of raspberry or ferns may reduce browsing compared to microsites dominated by mosses [29]. In heterogeneous stands, Rubus and other ground vegetation may vary in size and mask any effects of its quantity or quality on the level of browsing [20,21,23]. The height of the ground vegetation may partially explain the higher PB of medium-sized seedlings in comparison to small ones in our study, as the former standout of the vegetation for the whole year.
The connection between the presence of tree seedlings and browsing seems even less obvious. At low seedling cover, browsing tended to decrease with increasing cover, but at higher seedling cover, almost no effect or even the opposite was observed. In country-wide forest inventories in Austria and Switzerland, in which most plots contain only very few seedlings, the more browsing was found, the fewer seedlings were present [75,76]. In contrast, the incidence of conifer leader browsing by roe deer was positively associated with seedling density in plantation forests in England [55]. The relationship between browsing and seedling density may therefore indeed be curvilinear with much browsing at sites with very little regeneration, but also at sites with plenty of regeneration that stand out of the vegetation.

Conclusions
Quantifying ungulate densities via the numbers harvested and found dead by other causes within a certain forest area and period is a simple method frequently used in wildlife management in Switzerland. This approximation turned out to be highly positively correlated with the browsing probability (PB) for six frequent tree species in Switzerland. The ungulate density index (UDI, sensu [31]) that weighs different ungulate species had slightly higher explanatory power than single species densities. Moreover, sites in which some red deer and chamois were present had higher browsing on leader shoots than sites comprising only roe deer. Since the abundance of ungulates was one of the main factors influencing PB, focused hunting seems effective to reduce browsing in areas where tree regeneration is needed, in particular if applied over longer time spans, i.e., decades rather than just years [77].
However, browsing is a multifactorial phenomenon and ungulate density estimates alone only partly explain the variation in browsing probability. There are many other factors that also influence PB in tree seedlings of different tree species, including the important proximate driver's height class (with medium-sized seedlings more frequently browsed) and year. Also, distal contextual drivers such as the developmental stage of the stand, its basal area, vegetation cover, seedling cover and perhaps light availability and coarse woody debris affect the incidence of leader shoot browsing. Tree seedlings in small timber stands with intermediate basal area, a canopy openness of >20%, low overall vegetation cover, medium seedling cover and low amounts of fresh woody debris were found to be less browsed on their leader shoot. Although we included a large variety of variables that may influence the probability of leader shoot browsing, we neglected some potential drivers, including all landscape and land use factors (such as fragmentation, juxtaposition of feeding and cover resources, distance to pathways and roads, housing or human population densities). These factors could further explain the differences between the 49 indicator areas.
In conclusion, this study supports and summarizes-based on an extensive long-term dataset-previous findings found in many singular studies. As the factors driving browse risk are related to both forest structure and ungulate impact, forest and wildlife management are both important and should be coordinated in order to reduce browsing impacts.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/7/764/s1. Figure S1: Construction of canopy opening covariates, Figure S2: Distribution of the ungulate variables, Figure S3: Partial dependence plots with marginal effects of height class and year on browsing probability, Figure S4: Partial dependence plots with marginal effects of location and climate, Figure S5: Partial dependence plots with marginal effects of forest characteristics, Figure S6: Partial dependence plots with marginal effects of variables regarding ungulates, Figure S7: Partial dependence plots with marginal effects of location and light, Figure S8: Partial dependence plots with marginal effects of basal area and seedling cover, Figure S9: Partial dependence plots with marginal effects of vegetation cover, Figure S10: Variable importance in random forest with all plots of the 49 indicator areas, Figure  S11: Variable importance in weighted random forests, Figure S12: Partial dependence plots in weighted random forests, Figure S13: Further partial dependence plots in weighted random forests, Figure S14: Joint marginal effects of UDI and height class on PB, Figure S15: Subset analyses for ungulate composition.       Fixed effects are presented as odds ratio with 95% confidence intervals (95% CIs), random effects with the number of strata (N) and the variances.