Determining Statistically Robust Changes in Ungulate Browsing Pressure as a Basis for Adaptive Wildlife Management

: Ungulate browsing has a major impact on the composition and structure of forests. Repeatedly conducted, large-scale regeneration inventories can monitor the extent of browsing pressure and its impacts on forest regeneration development. Based on the respective results, the necessity and extent of wildlife management activities such as hunting, fencing, etc., can be identiﬁed at a landscape scale. However, such inventories have rarely been integrated into wildlife management decision making. In this article, we evaluate a regeneration inventory method which was carried out in the Bavarian Forest National Park between 2007 and 2018. We predict the browsing impact by calculating browsing probabilities using a logistic mixed effect model. To provide wildlife managers with feedback on their activities, we developed a test which can assess signiﬁcant changes in browsing probability between different inventory periods. To ﬁnd the minimum observable browsing probability change, we simulated ungulate browsing based on the data of a potential browsing indicator species ( Sorbus aucuparia ) in the National Park. Sorbus aucuparia is evenly distributed, commonly found, selectively browsed and meets the ecosystem development objectives in our study area. We were able to verify a browsing probability change down to ±5 percentage points with a sample size of about 1,000 observations per inventory run. In view of the size of the National Park and the annual ﬂuctuations in browsing pressure, this estimation accuracy seems sufﬁcient. In seeking the maximal cost-efﬁciency, we were able to reduce this sample size in a sensitivity analysis by about two thirds without severe loss of information for wildlife management. Based on our ﬁndings, the presented inventory method combined with our evaluation tool has the potential to be a robust and efﬁcient instrument to assess the impact of herbivores that are in the National Park and other regions.


Introduction
One of the most important processes in terrestrial ecosystems is the interaction between herbivores and vegetation [1]. Herbivores can influence vegetation at the individual, population and landscape scales [2][3][4][5][6]. In particular, large herbivore species such as ungu-change of browsing and, consequently, the expected mortality of the seedlings. In view of these correlations, Eiberle and Nigg [65] developed thresholds for permitted browsing intensity for common tree species in Swiss mountain forests. In general, these thresholds correspond to the previous average browsing intensity of a species within a height class in ungulate reach and the time that the seedlings take to cross this height class [17,60,65,66]. If the annually determined browsing intensity is exceeded on a long-term basis, the survival of this species is no longer assured. Shifts in the tree species composition in the subsequent stand can then be expected. Therefore, a time series analysis on browsing intensity, in alignment with browsing thresholds, indicates whether current wildlife management strategies have been successful according to set objectives (control function). Furthermore, they allow forecasting of the browsing development in the near future (planning function) [67,68]. Moreover, the continuous flow of information sensitizes the adaptive management about their impact, which mitigates the uncertainty of their future actions [48][49][50]69]. Nonetheless, according to Eibele and Nigg, "An undifferentiated assessment of browsing damage without checking its accuracy [. . . ] can easily lead to wrong conclusions" [65] (p. 776). The importance of this statement becomes especially clear in conditions with very low stem numbers or, rather, sample size. The browsing intensity does not provide confidence intervals to validate the quality and the descriptive power of the browsing information.
In summary, large-scale forest regeneration inventories have the potential to draw statistically robust landscape-scaled results of ungulate impact. National forest inventories, such as those described by Polley [70] or McWilliams et al. [71], take forest regeneration into account. However, either their temporal and/or spatial resolution are insufficient to allow wildlife management decisions to be made at a landscape or enterprise scale. With the exception of the inventory presented in this study, there is, to the best of our knowledge, no large-scale repeatedly executed regeneration inventory for monitoring landscape-scaled or enterprise-scaled ungulate impact. There also does not appear to be an instrument for displaying the inventory's accuracy.
The objective of this study is to identify a method for drawing statistically sound landscape-scaled wildlife management decisions based on large-scale regeneration inventories. The main focus will be on how inventory results can be placed into context with previous inventories to allow wildlife managers to adapt. We develop our statistical method in the context of the large-scale regeneration inventory applied in the Bavarian Forest National Park. Moreover, to be accepted and trusted by as many stakeholders as possible, we evaluate the inventory method applied in our study area. For the evaluation, we want to consider multiple aspects: First, the regeneration inventory method must be able to fulfill the inventory's purpose. Second, the implementation effort should be cost efficient [56,72]. A key factor in achieving a high level of cost-efficiency is having an optimal sample size to meet the inventory's purpose [73][74][75][76][77]. Third, the inventory design and data structure must be conceived in a comprehensible manner so that its results are independently reproducible to ensure comparability with other assessments [56]. It must also be repeatable so that trends can be identified through follow-up assessments to guarantee success controls [56]. To achieve our goals, we want to address the following questions: (i) How can the variability of browsing impact and its change over time be predicted and tested for statistical significance? (ii) Can the sample size of the established regeneration inventory be reduced and still serve wildlife managers as a valuable information source? (iii) What is the minimum observable significant change with respect to browsing impact in the Bavarian Forest National Park in applying its current forest inventory method?

Study Area
The data used in this analysis is based on a long term forest regeneration inventory in the Bavarian Forest National Park (BFNP). The national park covers 24,222 ha and connects with the adjacent Šumava National Park (68,064 ha) in Czech Republic to form the largest cohesive forest area in Central Europe-the Bohemian Forest Ecosystem [37,78]. Characterized by mountainous temperate forest ecosystems with Norway spruce (Picea abies) and European beech (Fagus sylvatica) as the most common tree species, the national park (NP) was established in 1970 and extended in 1997 [37]. The delayed amendment divides the BFNP into two large administrative zones: The "Rachel-Lusen" region (RLG), the initial area at establishment of the BFNP and the 10,903 ha large "Falkenstein-Rachel" region (FRG), which was used for forestry until 1997.
Despite anthropogenic influences on BFNP forests, forest communities are zoned by interactions of various site factors, such as height above sea level, climate and soil. Higher altitudes (above 1050 m) have shaped natural Mountain Spruce Forest Communities with over 90% Norway spruce. At the hillsides between 700 and 1250 m, the milder climate allows for mixed mountainous forest stands that are dominated by spruce (58%), beech (35%) and silver fir (Abies alba) (2%). Sycamore (Acer pseudoplatanus) and Norway maple (Acer platanoides), witch elm (Ulmus glabra), winter (Tilia cordata) and summer lime (Tilia platyphyllos), rowan (Prunus avium) and yew (Taxus baccata) are common secondary tree species. In valley areas (600-800 m) where cold and humid air accumulates and soils are waterlogged, the potential forest communities are limited to alluvial spruce forests. Norway spruce (83%) dominates the forest stands, followed by silver fir, beech and different birch species (Betula pendula and B. pubescens) [52,[78][79][80][81][82][83]. The areas of the high altitudes and valleys each account for 16% of the BFNP area and the hillsides for 68%.
While the proportion of spruce in the high altitudes and valley locations is quite natural due to the local site conditions, it was artificially increased over centuries by humankind, especially on the hillsides [78,81,83]. In particular, the share of silver fir, which naturally comprises one third of the forest in this ecological altitude level, was reduced [81]. The reduced diversity of forest composition, three storm events (1999, 2004 and 2007) and lack of human intervention due to the NP conventions incentivized a large-scale biotic ecosystem disturbance. In the two major outbreaks (1996-2000 and 2005-2009) of European spruce bark beetle (Ips typographus), over 6500 ha of spruce forest was infested [37,84]. It was shown that the die back of the mature spruce stands in the BFNP by bark beetle resulted in a fast onset of natural regeneration in areas of high density-especially the natural regeneration of the predominant spruce [85,86].

Wildlife Management
It is not only the forests but also the fauna of the BFNP that has been significantly changed by humankind until the establishment of the NP. Bison (Bison bonasus), moose (Alces alces), brown bears (Ursus arctos), lynxes (Lynx lynx) and wolves (Canis lupus) were exterminated, while roe deer (Capreolus capreolus), red deer (Cervus elaphus) and wild boar (Sus scrofa) populations increased through gamekeeping, improved living conditions and the lack of predation [87,88]. Since the 1990s, the lynx resettled in the BFNP from the Czech Republic and supported the regulation of ungulates, in particular, roe deer. In 2016, there was also evidence that a pair of wolves has returned to the BFNP [47,89]. Nevertheless, there is no evidence that the populations of red deer and wild boar are regulated by predators [52].
As the ungulates migrate in the winter months to areas with more favorable conditions (for example, less snow) which are partly outside of the BFNP's boundaries in managed forests, the BFNP administration is forced to regulate the ungulate populations within the NP to prevent browsing damage in their winter retreats [6,47,90]. The BFNP's wildlife management focuses on red deer and, in the past, also focused on roe deer [47]. Roe deer hunting was discontinued in RLG in 2007 to monitor roe deer population dynamics and to leave the deer population to its natural predation pressure. In 2012, roe deer hunting was also halted in FRG. The long-term inventory of ungulate browsing on the forest regeneration is an integral part of the NP's adaptive wildlife management. Conclusions about the state of regeneration and the development of the ungulate densities can be drawn from the browsing changes over time [52].

Inventory Method
Annual monitoring of the regeneration and the ungulate browsing inventories has taken place in RLG since 2007. Since 2008, the FRG was also included in the regeneration inventory. In 2012, the inventory frequency was reduced to a 3 year interval.
The inventory method, including sampling design, is based on the Bavarian regeneration inventory [91][92][93]. The Bavarian regeneration inventory is designed to monitor the seedling state on forest regeneration sites so that silvicultural regeneration objectives can be achieved [91,92]. Only regeneration sites within forests (natural regeneration as well as artificially established seedlings) are recorded during the inventory process.
An 800 × 800 m (64 ha) mesh grid spanned across the entire BFNP area (see Figures 1  and A9). Starting from the grid point, the closest regeneration area is identified. To avoid double counting, the regeneration area size should be at least 50 m in diameter and should have a regeneration density of more than 1300 plants ha −1 ≥ 20 cm. It was ensured that the leading shoots of the observed seedlings are still within reach of the locally occurring herbivores. In the BFNP, the maximal possible browsing height (considering the snow conditions in winter) is set at 200 cm. Measurements only take place on regeneration areas that meet these criteria. There is no grid where, due to the aforementioned requirements, no regeneration information is available. To cover spatial browsing variation inside a regeneration area, a transect of at least 40 m (maximum 100 m) length is laid through the widest section of the regeneration area. Five sample points are mapped on the straight transect at intervals each corresponding to a quarter of the transect length. The first and the last sample points have to be at least 5 m away from the edge of the regeneration area. The coordinates in the center of the first sample point are recorded. Together, all sample points of a regeneration area form a (cluster) plot (see Figure 1). Each cluster plot is individualized by a plot ID [52,93]. Starting from the center of each sample point, the 15 closest individual seedlings are examined, which fulfill the height requirements (20-200 cm). The distance of the 15th recorded regeneration plant to the center of the sample point corresponds to the sample points' radius. On each seedling inventoried, tree species and height are determined. To monitor the browsing injuries on the terminal shoot that occurred since the beginning of the last vegetation period, binary notation was used.
This standard procedure was extended by a separated regeneration inventory to increase the number of samples of rare tree species [94]. Tree species are classified as "rare" if they account for less than 5% of the total regeneration. In the BFNP, these are all of the tree species present with the exception of the Norway spruce and European beech species. At the centers of the sample points 2 and 4 of each transect, up to five of the closest "rare seedlings", which are ≥20 cm and smaller than the maximum browsing height, are mapped at a maximum possible sample radius of 5 m.
Until 2010, about 400 plots per inventory were recorded. However, under the consideration of the inventory extension for rare tree species, the plot quantity was reduced to ca. 250 from 2011 onward [95].

Data Set
The data set contains the year of the respective inventory, the section of the BFNP, plot ID, consecutive number of sample points (1)(2)(3)(4)(5), consecutive number of observations (1-15), tree species, seedlings' height and terminal shoot browsing (TRUE/FALSE).
We merged the data set of the standard regeneration inventories with the inventory extension. We took care to exclude double counting of rare tree species at the plot's sample points 2 and 4. Therefore, we removed the rare tree species (all except spruce and beech) at the above-mentioned sample points from the standard inventory in all years from 2010 onward before merging the extended data. In total, there are 183,452 observations in our data set.
We took all tree species with a relative proportion of more than 0.5% in the data set (incorporating the inventory extension) into account in our analysis (compare, Appendix A.1). The tree species range includes Norway spruce, European beech, rowan, silver fir, birch (Betula spp.) and sycamore maple. Below this threshold, the quantity of seedlings was not sufficient to reliably determine browsing changes over a longer period of time.

Browsing Probability and Its Change
Ungulate impact was predicted by estimating the odd-ratios of a single tree being browsed by using a generalized logistic mixed effect regression model (GLMM)-the browsing probability (BP) (question i) [6,26,96]. The logistic model used was chosen due to the binary nature of the single-tree browsing observations and to ensure that the resulting model estimates values between 0 and 1. We decided to apply a mixed model since the browsing observations clustered by every plot over time are expected to be stochastically dependent. As a fixed effect, we included year (of inventory) as a categorical variable (factor) with a fixed intercept and set a random intercept for the plot ID (see (1). The random intercept enhances the logistic model by adding a measurement repetition factor [97][98][99]. Furthermore, the random intercept serves as a surrogate for the possible effects of covariates that were not considered during data acquisition and, thus, not measured [97,98,100]. For example, browsing might depend on the altitude level and a plot can only be assigned to exactly one ecological altitude level. Such unobserved heterogeneity is therefore taken into account in our GLMM. In a previous model, we tried to account for possible BP variability between sample points nested in a plot in the random intercept. However, since this resulted in a decrease in model quality (Akaike information criterion) and since no or only a little variability was shown between sample points for many tree species, we have not pursued random intercept modification any further.
where η ij = linear predictor for the browsing probability; β 0 = fixed intercept (the origin in our model); u i = cluster-specific (random) deviation from the fixed intercept; β 0 + u i = random intercept for cluster (plot) i; β 1 = fixed slope parameter of covariate t j (year); ij = error term; P(y ij = 1) = browsing probability.
We estimated the BP for every tree species with the data set of the whole area of the BFNP and with the data set divided into its two management sections (FRG and RLG). Plot-level estimates are also potentially possible. Nonetheless, our analyses are carried out for the landscape and management unit scale since wildlife management decision making at the plot level would be impractical.
Subsequently, we formulated linear hypotheses for a simultaneous ANOVA procedure for each regression coefficient of the year j (see Equation (3)), tested them for significance and computed their 95% confidence intervals [96,101].
To estimate the changes in the BP based on the GLMM (question i), we modified the linear hypotheses (see Equation (4)). For each hypothesis, we took the difference between the regression parameters of the year j and the previous year of inventory j − 1. Testing the BP change over other time periods is, of course, also possible.
However, since the regression parameters are subtracted as logits, converting the BP change into a response (percentage) value afterwards is not possible.

Sensitivity Analysis
In seeking cost efficiency for the inventory method (question ii), we first determined the least number of observations of a tree species that must be present in the data set to maintain a significant change of BP between two years. This requires that a change was classified as significant in previous calculations. First, we divided the data set into two years in which the significant change was observed. Then, we ran a sensitivity analysis by using a Monte Carlo simulation, which iteratively reduced the total number of all seedlings across tree species. To achieve this, we simultaneously removed 20 random sample points in each iteration step from the data set between the years to be compared. Afterwards, we estimated the BP change with its confidence intervals for the respective tree species. We repeated this operation 1000 times for each step, to ensure that the estimations were based on different sample point combinations. For evaluation and graphical representation of the simulated values, we took the means of the logarithmic BP change in each iteration step (and their confidence intervals). The resulting minimum sample size marks the point at which the browsing change was last classified as significant. Between two data acquisition periods, the quantity of observed seedlings usually fluctuates. This difference between two years cannot be compensated by the simultaneous reduction in the samples. The minimum sample size always corresponds to the smallest sample size between each year.
In the results and appendix, only the tree species that represent more than 0.5% of the data set and that possess a significant BP change between the two most recent pairs of inventory years can be found (2012-2015 and 2015-2018). Assuming that the browsing probability estimates are more accurate due to a larger sample size, we used the data set of the BFNP as a whole.

Simulation of Ungulate Browsing
In the sensitivity analysis, we only worked with the assessed inventory data. Finding the minimum sample size for a significant BP change of a tree species was, therefore, limited to the actual observed browsing change. Theoretically, it could be possible that the applied inventory method is capable of detecting even smaller BP changes than observed in nature. To address question iii-to determine whether, for example, an increase in BP of exactly 2 percentage points (PP) for a tree species would represent a significant change-we generated data sets by simulating ungulate browsing that exactly reflect a certain change.
For this simulation we selected rowan because of its even distribution in the NP (throughout all altitude levels), its high amplitudes of BP and because it has the largest sample size among the rare tree species (compare Sections 3.1 and 4.4). We used the most recent data from 2015 and 2018 for the simulation, as it represents the latest spatial distribution and abundance of rowan.
We left the browsing impact of the 2015 data set unchanged. For the year 2018, however, we set all browsing observations to "FALSE" to obtain a BP of 0%. In each simulation step, the simulation process randomly selects from the rowan 2018 data pool of five individuals. The selected individuals are then "browsed" by assigning the logic value "TRUE". We then estimated the BP. As soon as the BP increased by ≥1PP compared to 2015, the data set, with the individuals browsed by the simulation, was extracted and saved. This simulation sequence repeated for a BP increase of +1 to +15 PP in single absolute steps. The same sequence was modified to simulate a decrease in BP of −1 to −15 PP. Afterwards, we checked the change of BP for significance based on the generated data. The smallest and still significant change then equals the estimation accuracy for rowan of the currently applied inventory method in the BFNP. Finally, we applied the sensitivity analysis of Section 2.5.2 to the extracted data sets which were tested as significant. We could then examine the minimum number of observations needed for different estimation accuracies.

Computational Details
All calculations and simulations were executed by using R version 3.5.3 [102]. For the GLMM executions, we used the "lme4" package [98]. For testing the multiple linear hypothesis in consideration of simultaneous interference, we used the "multcomp" package [96]. All necessary editing of data structures and the graphical evaluation of all results of this work were performed with the package "tidyverse" [103].

Browsing Probability Time Series
The detailed results of the BP and the BP change of the logistic mixed model of the respective tree species and NP region can be found in the Appendix A. To focus this work on the evaluation tool, only the highlights of the BFNP's browsing probability times series are presented.
The BP levels between data acquisition years and tree species differ substantially. While the BP of the most frequent tree species (Norway spruce) never exceeded 1.1% (see Figure 2A), much higher BP levels and also BP amplitudes were found for rare broad-leafed tree species (see Appendix A). For instance the BPs of rowan berry fluctuate between 10.5% (2009) and 34.9% (2015) (see Figure 2A). Furthermore, we noted a trend that the rarer a tree species occurs in the NP (or the smaller the sample size), the less accurate the BP estimates become (or the wider the confidence intervals) (see Figure A4A). With few exceptions in tree species with a sparse sample size, the categorical expressions of variable year had a significant effect on the BP estimates (see Table A1).
To raise certainty that a BP change happened on the population level, we tested the logarithmic change of BP for significance (see Figures 2B and A4). We found that the larger the sample size and the greater the BP change between years, the more likely a significant BP change could be detected (see Appendix A and Table A2). For instance, the lowest observed but still significant BP change in the entire NP is for rowan, which had an increase of +4.8 PP. For silver fir, although its sample size is slightly smaller than rowan, a BP change of +4.2% was observed (see Figure 2B). However, silver fir (ca. 8.5 obs. per plot with fir in 2018) is spatially more clumped compared to rowan (ca. 8 obs. per plot with rowan in 2018), which could explain the described phenomenon (compare Figure A10).

Significance of Changes and Sample Sizes
The simultaneous sample size reduction process demonstrated its flawless functionality. The smaller a sample size of a tree species became between two years, the more likely the once significant BP change became insignificant. We can reemphasize that the larger the sample size and the larger the observed significant BP change, the more observations could be removed from the sample size (see Appendix A.2). For example, we were able to reduce the sample size of Norway spruce by over 50% of its recorded sample size without losing significance for a BP change of +0.5 PP between 2015 and 2018 (see Figure A8). For the rare rowan berry tree species, we were able to reduce the sample size by ca. 80% for a BP change of +19.5 PP between 2012 and 2015 (see Figure A7).

Simulation of Ungulate Browsing on Rowan
We applied the browsing simulation on rowan only. Rowan is the most common rare tree species in the BFNP that sensitively reflects changes in ungulate browsing and it is spatially evenly distributed with a sufficient sample size (see Sections 2.5.3 and 4.4.1 and Figure A10).
By increasing step-wise the BP of rowan in the 2018 data set by one percentage point each, starting from 34.9% (2015), we found a significant change starting from an increase of +5 PP (see Figure 3). The same BP change also results in a significant difference when we lowered the BP. We saved and provided data sets with the increased or decreased sample size for the following sensitivity analysis. The data sets contained the significant changes between ±5 PP and ±15 PP.
We were able to reduce the sample size for a significant increase in BP of +5 PP from 1007 to a minimum of 757 observations. A significant BP decrease of −5 PP required at least 882 observations. A minimum sample size of 399 (increase) or 382 (decrease) was required to observe a BP change by ±8 PP. For a significant BP change of ±15 PP, 192 observations (increase) or 191 (decrease) were required (see Figure 3).
The initially rapidly decreasing minimum sample size with increasing BP change must be emphasized. From a change of about ±9 PP, the decrease in the minimum sample size mitigates and remains relatively constant at a level above 200 observations. We were able to fit an exponential function from the minimum required sample sizes per BP change with a coefficient of determination of over 90% (see Figure 3). The estimated graphs of Figure 3 reflect their inflection point at a minimum sample size of ca. 350 observations. The exponential function for a significant BP decrease is listed in Equation (5) and for a BP increase in Equation (6).

Evaluating Ungulate Browsing Impact (Question I)
The results show that calculating the BP with confidence intervals instead of the browsing intensity provides adaptive wildlife management with a better understanding of the accuracy of browsing predictions. Similar but rare approaches can be found by Kupferschmid et al. [26] and Hothorn et al. [96]. Morellet et al. [54] or Cederlund et al. [104], however, pointed out the importance of observing ecological changes for managing herbivore populations. A population level BP change, though, cannot be detected with certainty just by confidence interval widths of estimated BPs. Certainty about change is crucial when evaluating the success of wildlife management activities, e.g., with browsing thresholds [49,65,69].
Our example data set consists of 65% spruce observations. The estimates of BP and BP changes are correspondingly precise. Considering the low spruce BP compared to their browsing threshold of 12% [65], the estimation accuracy of the BP change is disproportionate to the effect it has on the spruce population. However, the closer a tree species reaches its respective browsing threshold, the more important it becomes to determine with certainty whether it exceeds it. For instance, we can observe in Figure 2A that the confidence intervals for the BP predictions of silver fir for the years 2011 and 2012 overlap. However, the significance test ( Figure 2B) of the logarithmic change detects a significant reduction. While the 2011 BP estimate for silver fir exceeds the fir's critical browsing threshold (of 9%), the 2012 estimate is below it. The silver fir example underlines the importance of a differentiated evaluation of statistically significant ungulate browsing changes to prevent false conclusions, enabling adaptive wildlife management.

The Ideal Sample Size and Cost Efficiency (Question II and III)
An ideal sample size depends on the resources of the organization conducting the inventory [105,106], the spatial distribution of the object of observation [107,108], the BP level [109] and the inventory objective. The sample size possesses indispensable information for the purposeful planning of sample-based inventory procedures [77], as it is decisive for the accuracy of an estimate. We measured the estimation accuracy in our results by the possibility of identifying a change of BP as significant. Our results outline that a higher sample size does not generally equal a higher estimation accuracy. For the spatially clumped silver fir, which has a slightly lower sample size than rowan and considerably lower BPs, we observed a significant BP change earlier than for rowan (see Sections 3.3 and Appendix A.2.2, Figure A10).
As the estimation accuracy is altered by the site conditions in the inventoried area, the inventory method has to be designed so that an accuracy can be reached, which is suitable for the inventory's objective-in our example, the state of forest regeneration [93]. The estimates must be precise enough to match browsing predictions and changes with the browsing thresholds of the tree species in focus [54]. Furthermore, a sufficient degree of estimation accuracy, especially for more rare tree species, is commonly counteracted by the limited resources of the inventory taker [105,106]. Therefore, inventory planning should always aim to achieve the maximum possible cost efficiency [105]. This can be accomplished either by achieving the highest possible accuracy with a fixed budget or by reducing the effort to the lowest possible level with fixed accuracy (compare [105,106,110,111]). Accordingly, we strive to find an ideal minimum sample size (as a guideline) for a desired level of accuracy based on the natural distribution of a species in our study site the BFNP.
Considering that the selectively browsed rowan could be a suitable browsing indicator tree species reflecting the management objectives of the BFNP (see Section 4.4.1), we based our browsing simulations and sensitivity analysis on its data set [68,112]. The results should provide an example of application of the simulation tool by using the most recent data from 2015 to 2018. We found the most cost-efficient minimum sample size of rowan to be ca. 350 observations for the BFNP with an accuracy for detecting changes between [−9:−8] and [8:9] (see Figure 3): Within this range, the minimum sample size decreases most effectively compared to the estimation accuracy since it marks the graph's inflection point. In view of the annual BP fluctuations of rowan, the estimation accuracy of this cost-efficient sample size still seems sufficient.
However, as the ideal sample size is influenced by several conditions depending on the inventoried area, we recommend defining an inventory objective to which a desired estimation accuracy is linked. In the long term, the sample size could then be optimized to reach cost-efficiency and to define an ideal sample size with a comparable sensitivity analysis.

Inventory Design and Further Improvements
The data acquisition method used in the context of our study is based on the widely applied point-to-tree distance sampling technique (also known as k-tree or fixed count sampling) by Hothorn and Müller [6], Prodan [113], Kleinn and Vilčko [114], Staupendahl [115], Lynch and Rusydi [116], Picard et al. [117]. Distance sampling is known for its ease of application in the field; the amount of work remains constant at high stand densities, which are common on regeneration sites [115][116][117].
However, when working on distance sampled data, statistical challenges arise. Unbiased estimates of stand or regeneration densities require that the selection probability of all surveyed population elements in a plot design to be known [114,118]. The selection probability of an individual, using a distance sampling technique, depends on knowing the spatial distribution of its surrounding individuals [114,119]. Efforts required for calculating density-dependent selection probabilities of every acquired individual would make distance techniques impractical and are not usually provided [114,119]. Without incorporating selection-probabilities, absolute stand densities will be overestimated [114]. Relative browsing impact estimates are unlikely to be affected by sampling design bias. Nevertheless, the larger k is, the smaller the bias of estimates based on distance sampled data [114]. Therefore, a bias caused by distance sampling with k = 15, as in our case, is presumably not relevant and is outweighed by the ease of sampling.
However, the inventory method applied excludes all regeneration sites with a density lower than 1300 seedlings ha −1 . This sets the selection probability of individuals occurring in areas with low regeneration density to zero. Kuijper et al. [120] and Churski et al. [121] both observed that ungulate browsing is clustered on high density regeneration sites in forest gaps. Consequently, the ungulate impact on the total regeneration of the BFNP might be overestimated with the current regeneration site selection rule. At the same time, the seedling selection is limited to individuals above 20 cm. The browsing pressure and the shift of competitive conditions of selectively browsed species starts from a seedling's germination [6,[122][123][124][125]. Therefore, the BP estimates could be left truncated. Moreover, the variable transect size (which varies between 40-100 m) alters the selection probability of seedlings again. Furthermore, the inventory method excludes the edges of regeneration sites due to its 5 m wide buffer zone (compare [93]). Clustered browsing on the edge of regeneration sites with a high density, such as observed byČermák et al. [126], cannot be captured. The selection probabilities of rare tree species are further increased by the inventory extension. Hence the selection probability of a seedling, which is already hard to grasp when using a distance sampling technique, is further modified by multiple additional constraints which have only minimal practical advantages. In addition, the inventory extension falsifies statements about the relative species composition by increasing the sample size of rare tree species. Therefore, we could not consider the selection probability for the purposes of the application of our evaluation tool, for estimating BPs or their change applied in the BFNP example.
It is likely that our relative BP predictions are not biased by the applied distance sampling technique itself. Nevertheless, the inventory design is limited by the summed constraints. These constraints result in an unpredictable bias, reduces the comprehensibility of the inventory method and, therefore, reduce its acceptance. These issues could be easily resolved by fewer or no constraints on the inventory method. We want to highlight that the inventory method is only suitable to monitor ungulate browsing pressure; unbiased regeneration densities cannot be predicted, it only provides an impression about the proportions of the species composition and wildlife management must interpret the BPs as the upper limit of the observed area, as browsing pressure is likely to be higher on denser regeneration areas. Despite this criticism, regeneration inventories supporting landscape-scaled or enterprise-scaled wildlife management remain a rarity. In most cases, wildlife managers cannot refer to any data basis. Even if, from a statistical perspective, the inventory method requires improvement, it was sufficient from a practical stand-point for evaluating the ungulate impact on denser forest regeneration sites, thus, supporting adaptive wildlife management. The results of the BP time series evaluation show how diversely individual tree species are affected by browsing (see Section 3.1 and Appendix A). Therefore, it is advisable that wildlife management decisions are based on browsing information of tree species that meet the overall management objectives of the considered forest area. These objectives may, for example, reflect forest owner-specific management interests [35,127] or be legally required, as in the BFNP. The most precise BP estimates in the BFNP example were obtained from spruce, which accounts for 65% of the observations in 2018. Due to its low palatability [19], spruce is not selectively browsed and, therefore, does not reflect the current browsing situation in the NP in terms of the total regeneration. For an inventory that aims for cost-efficiency, the high data acquisition effort cannot be justified for the evaluation of the browsing impact on spruce. Sycamore maple, on the other hand, is a selectively browsed tree species [125,128] and cannot provide reliable browsing results in the NP due to its low abundance of less than 1% in the data set. Alternatively, instead of increasing the global sample size, species could be identified that are selectively browsed and meet the inventory objective. In other words, browsing indicator species could be determined.
Even though an indicator species analysis is a common feature of ecological research (e.g., [129,130]), we were able to find just a few examples in forestry literature within the browsing context. For instance, Anderson [131] used Trillium grandiflorum as an indicator species to show the browsing intensity of the white-tailed deer (Odocoileus virginianus). Chevrier et al. [51] have observed how browsing intensity in young selectively browsed oak species (Quercus spp.) is related to the roe deer density in a 10 year series of experiments in France. From their observations they established an "oak browsing index", which increased linearly with a high correlation to increasing roe deer densities. Furthermore, Morellet et al. [54] identified several indicator species to infer the size and performance of roe deer populations and their impact on the habitat to support an adaptive wildlife management.
However, using indicator species may result in a bias in favor of them, at the expense of other under-represented species [132]. To counteract this phenomenon, multiple requirements must be met to generate valuable information for adaptive wildlife management from indicator species. The following characteristics, which are modified and applied to wildlife management purposes according to Anderson [131], Carignan and Villard [106] and Williams et al. [133], stand out. The indicator species should: • Act as an "early warning system" that reflects the influence of ungulate browsing of other tree species or, rather, a development of the ungulate densities in time; • Be preferred by ungulates so that the indicator species sensitively indicates changes (large amplitude) rather than just indicating a change. Ideally, the species should be so sensitive that strong browsing quickly reduces its abundance; • Be able to be monitored cost-efficiently; due to a high abundance and uniform distribution in the study area, an ideal browsing indicator species achieves a high estimation accuracy in data collection without great effort.
For the BFNP, rowan could be used as a well-fitting indicator species. Not only is a sufficient number of rowan available in the NP, but this species is also distributed evenly throughout the park to representatively observe browsing changes (see Figure A10). Rowan also seems to show changes of the browsing impact in the BFNP sensitively, which can be observed on its broad browsing amplitude (see Appendix A.1). Furthermore, rowan is considered a pioneer tree species and a biodiversity carrier [68,134] and, thus, it could meet the management objectives of the BFNP [46].
However, further application studies are needed to investigate the effect of adaptive wildlife management based on an indicator tree species on wildlife populations and forest regeneration. Moreover, inventory methods could be sought to efficiently study a indicator tree species.

From Which Browsing Probability or Threshold Does the Population of a Tree Species Decrease in the Long Run?
The correlation between browsing, loss of height increment and mortality has already resulted in the development of thresholds for the long-term permissible browsing intensity. These thresholds serve as a reference value for the respective tree species for wildlife management and are, therefore, of major importance (e.g., in the BFNP). However, the thresholds of Eiberle and Nigg [65,128] were determined based on 1363 observations in total, recorded at 23 different sites of the montane vegetation zone in the Swiss Alps. Accordingly, the estimated thresholds are only valid for this habitat and are not universally applicable without restriction. Additionally, they are limited to a small tree species spectrum. To the best of our knowledge, there is no further research to extend or update these thresholds. Only Kupferschmid et al. [25] recently documented a significant reduction in seedling stem numbers from a browsing probability between 5-10%. These observations relate to a wide range of different tree species that are also grown in Swiss mountain forests. Without comprehensive reference values, a regeneration inventory generates data without being able to evaluate it.
Light is often the most limited resource on the forest ground [39,135]. Compared to light, soil nutrients as well as water play a less important role for the growth of seedlings and saplings in most cases and in later stages in the temperate zone [3,136]. Under low light conditions and when the seedlings are already stressed, it is often not the one-time browsing that results in death but the repeated removal of the terminal shoot [19,65,[137][138][139]. Nonselectively browsed seedlings gain an advantage in the competition for scarce light and can outcompete the browsed individuals-the so-called apparent competition [140]. Even if the primary shoot of a successively browsed tree grows out of ungulate reach, the competitive shift renders a subsequent death more likely [140,141]. These deer-plant, plant-environment interactions and tree-population dynamics should be taken into account when creating new browsing thresholds [60,64,135,142]. For this purpose, the development of sample trees, which forms the basis for the calculation of thresholds, has to be recorded spatially. The sample trees have to be placed into relation to the respective tree species and their height [64,135,142]. Such data acquisition, however, requires very high levels of effort since a very large spectrum of naturally growing species and mixture ratios has to be recorded over years. New three-dimensional forest development modeling approaches could help to better understand the tree population dynamics in combination with ungulate browsing (e.g., [135,[141][142][143][144][145]). Individual survival probabilities (the likelihood to grow to an adult height) could then be calculated for specific tree species, affected by certain browsing probabilities and species constellations on a regeneration site.

Conclusions
Large-scale regeneration inventories have rarely been implemented in local adaptive wildlife management concepts and have seldom been the subject of research up to now. The Bavarian regeneration inventory, which was with minor modifications executed in the BFNP, is one of the first of its kind. Although it has room for methodical improvement, after little methodical adjustments of the inventory procedure, we can recommend the applied regeneration inventory as a base for adaptive wildlife management.
Browsing probabilities (BP) and their change computed by a generalized mixed effect model from periodically acquired regeneration data provide valuable information for wildlife management and underline the necessity of such inventories. Furthermore, estimating BPs instead of calculating the browsing intensity can significantly increase the credibility of regeneration inventories. Confidence intervals, which are mostly affected by the sample size and the local variability of the browsing impact, can be calculated, which further allows the detection of significant changes over time. BP changes can then be compared to browsing thresholds to estimate the effect of the change on the population of the tree species in focus. Our browsing time series clearly shows how much and how rapidly the ungulate impact can fluctuate. It also demonstrates how estimation accuracies can differ between tree species. A high degree of validity and credibility of the regeneration inventories are essential for wildlife management and for its public acceptance. Ungulates, which are managed based on such inventory reports, are a natural component of forest ecosystems; the forest is their livelihood. A regeneration inventory should be sufficiently precise to find a balance between forest regeneration and ungulates, improving the living conditions of the animals without compromising the objectives of silviculture and/or nature conservation.
The sensitivity analyses and the browsing simulations applied on a potential browsing indicator species show, that sample size could be reduced in favor of cost-efficiency and are still able to retrieve robust statements about the current browsing situation and its development. The analyses not only proved its high potential as a tool for saving inventory costs in the long term, but also demonstrated that an ideal sample size depends on several site dependent factors. Therefore, instead of defining an ideal or minimum sample size for targeted species, which is variable in every use case, we suggest determining a minimum estimation accuracy. A minimum estimation accuracy combined with an unambiguously formulated inventory aim would significantly increase the information value of a regeneration inventory and serve adaptive wildlife management as a reliable resource.  Acknowledgments: We thank the program Ziel ETZ Free State of Bavaria-Czech Republic 2014-2020 (INTERREG V) and the program "Biologische Vielfalt" for funding the present research. We would also like to thank Sennhenn-Reulen, who advised us and provided valuable tips with statistically driven questions. We are grateful to all the participants who acquired the browsing data in the Bavarian Forest National Park. Finally, we thank Alena Anderson Chilian and Juanita I. Schmidhammer for language editing and the anonymous reviewers for their valuable suggestions.

Conflicts of Interest:
The authors declare no conflict of interest and the funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: Norway spruce is the most common tree species in the regeneration data set. In the most recent regeneration inventory from 2018, 10,774 individuals were recorded throughout the NP (see Figure A1). Between 2010 and 2011, we observed a strong decrease in spruce individuals in both NP areas in our data set. This change is due to the optimization and thus reduction in the plot number in the regeneration inventory (see Section 2.3). However, the relative shares of spruce in relation to the total amount of regeneration in the NP in the respective year decreased from about 69% (2008) to 58% (2018).   The BP over the entire NP corresponds to the average BP of RLG and FRG (see Figures A4A, A5A and A6A), although the RLG is slightly more dominant due to its sample size (see Figures A1 and A3). The time series starts with a 0.3% BP identical to RLG since no inventory data were available for FRG in 2007 and remained relatively constant from 2008 onwards at about 0.9%. The strong increase in browsing probability in 2015 in FRG was also reflected in the evaluations of the entire NP. In 2015, spruce browsing culminated at approximately 1.1%. Moreover, the confidence interval width increased to ±0.3 PP, while the average width of all years was ±0.2 PP. In the most recent year, spruce BP was 0.6% with a 95%-confidence interval width of ±0.2 PP. Compared to all other tree species, spruce had the lowest browsing probability.
To raise certainty that a BP change happened on the population level, we tested the logarithmic change of BP for significance ( Figure A4B). For spruce, all changes of the BP are significant apart from the difference between the years 2010 and 2011. The smallest significant change for the whole NP was +0.3 PP from 2009 to 2010.     With a share of 24%, beech is the second most common tree species in the regeneration data set (see Figure A1). The evaluations from 2018 were based on 5250 beech seedlings and seedlings of which about one-third were in FRG and two-thirds in RLG (see Figures A2 and A3). As in the case of spruce, a lower number of beech trees were observed from 2010 onward due to the reduction in the plot numbers (see Figure A1).
The estimated BP of beech is considerably higher compared to spruce. BP had its minimum of 3.5% in 2009 and its maximum of 13.3% in 2015. While the confidence intervals in 2015 ranged around the estimate by ±2.7 PP, they were ±1.5 PP wide on average (see Figure A4A). From 2010 to 2011 (+1.4 PP) and from 2011 to 2012 (±0.0 PP), we could not detect any significant browsing changes in the entire NP area. The lowest observed and still significant change of BP was found between 2009 and 2010 with an increase of +1.4 PP (see Figure A4B).

Appendix A.1.3. Rowan Berry
The rowan is the most frequently observed rare tree species over all years. Observations numbering 1251 were available in 2018, which make up 6.7% of the data set. Only in 2007 and 2015 was silver fir the more frequent rare tree species in the NP (see Figure A1). In 2018, the rowan had about the same share in the regeneration data set in both NP subregions (see Figures A5 and A6). The inventory extension for rare tree species from 2010 onward (see Section 2.5.1) and the optimization of the plot number are quantitatively noticeable in the rowan data set. Due to the inventory extension, the share of rare tree species in the total regeneration data is distorted so that we cannot observe any trends here.
Rowan was browsed significantly more than spruce or beech. Between 2007 and 2012, the BP fluctuated between 10.5% (2009) and 18.7% (see Figure A4A). From 2010 onward, we noticed a continuous increase in BP, culminating in 2015 with 34.9%. Here, the relatively low BP of 28.1% in RLG in 2015 could compensate for the high BP in FRG of approximately 47% (see Figures A5 and A6). In 2018, we estimated the probability of a rowan to be browsed at 33.9%. At about ±4.1 PP, the confidence intervals ranged around the estimates.
In With around 1062 observations in 2018, silver fir is the second most common rare tree species in the data set (inventory extension included (see Figure A1). The relative share of fir in the regeneration data set has constantly been rising from 2010. Nevertheless, it does not correspond to the share of actual regeneration since the inventory extension raises the fir sample size (see Section 2.5.1). In direct comparison, significantly more fir trees were mapped in RLG in 2018 than in FRG-823 (77.5%) to 239 (22.5%) (see Figures A2 and A3).
Up to and including 2012, the BP of fir in the entire NP remained below 10%. The BP more than doubled from 5.1% to 11.2% only from 2012 to 2015. Due to the increase in BP in RLG towards 2018 (11.5%) and the already high values in FRG (27.6%), a further increase to 15.7% was observed in the entire area from 2015 to 2018 (see Figures A5 and A6). On average, the confidence intervals surround the BP estimate ±2,6% (see Figure A4A).
From 2011 on, we were able to address all BP changes with a safe 95% confidence interval for the entire NP as well as for the RLG. The smallest and still significant BP change is +4.2% (from 2012 to 2015), while in the RLG, a change of +2.9% (from 2012 to 2015) is still significant despite a smaller sample size(see Figures A4B, A5 and A6).
Birch trees are the third most common rare tree species group (including the inventory extension) in the NP's data set. Observations numbering 241, which correspond to a share of 1.3% in our data set, are included for our BP estimation in 2018. Of the 241 seedlings (2018), 210 (87%) were in RLG and 31 (13%) in FRG (see Figures A1-A3). Due to the spatial dominance of birch in RLG, the BP estimates for the total NP largely correspond to those of the RLG. For the same reason, it is not possible to make any statements about the development of browsing in FRG.
Ungulate browsing on birch was strongly fluctuating over time in RLG. In 2007, the time series started with the maximum BP value of 32% and decreased to 11.1% in 2010 but doubled in 2011. In 2012, we observed a strong decrease in the browsing probability to 1.6% and an increase in 2018 to 26.3%. Similar to the fluctuating BP, the confidence interval widths were also fluctuating between ±3.3 PP (2012) and ±19, PP (2008) (see Figures A4A,  A5 and A6).
The determined logarithmic BP changes in the RLG and the entire NP show that only the changes from 2011 to 2012 (−19.2 PP) and the increase in ungulate browsing by +13.1 PP from 2012 to 2015 are significant (see Figures A4B, A5 and A6).
Appendix A.1.6. Sycamore Maple Sycamore maple is the fourth rarest tree species in the NP (with inventory extension). With 67 observations, the tree species accounts for about 0.4% of the examined regeneration in 2018. The number as well as the percentage of observed sycamore trees seemed to decrease until 2012. In contrast to the rare tree species described above, the share of sycamore maple in the regeneration data set of the respective sub-area is balanced at approximately 0.4% in RLG and FRG in 2018 (see Figures A1-A3).
Between 2008 and 2018, we estimated sycamore maple probabilities of BP in the NP area between 5% (2009) and 25.3% (2018). From 2012 onward, we identified a steadily increasing BP. On average, the confidence interval width was about ±12.3 PP wide (see Figure A4A). Only the increase in BP in the entire NP by +12 PP from 2009 to 2010 was significant. While we also measured the change between 2008 and 2009 in the RLG as significant, we did not observe any significant change in the FRG (see Figures A4B,  A5 and A6). We confirmed the increase in spruce's BP of +0.6 PP in Appendix A.1.1 as significant. The estimation was based on 12,073 seedling observations in 2012 and 12,280 in 2015. We were able to reduce the sample size down to 5008 individuals without losing significance once, which was observed at the excess of the zero threshold of the confidence interval boundaries. Below 5008 observations, we sometimes observed a loss of significance. From approximately less than 3600 individuals, the significance loss incidence becomes regular (see Figure A7).
Between 2012 (4680 observations) and 2015 (4615 observations), the BP of beech rose significantly by +8 PP (see Appendix A.1.2). The observed change was so distinct that the sample size reduction was not able to break the significance. If we reduced the sample size below 992 beech seedlings, the sensitivity analysis was terminated by multiple error messages, which implied that the sample size per plot was too low (see Figure A7).
Between 2012 (1315 observations) and 2015 (1007 observations), we observed a change of BP of rowan by +19.5 PP (Appendix A.1.3). Within the sensitivity analysis, we were able to reduce the minimum sample size to 198 without losing significance until the analysis process was aborted by the same error messages as mentioned above (see Figure A7).
For silver fir, BP doubled from 2012 (799 observations) and 2015 (1053 observations) by +6.1 PP (see Appendix A.1.4). We were able to observe this significant change until a minimum sample size of 189 seedlings (see Figure A7).
Between 2012 and 2015, the estimated BP of birch increased six-fold from 2.7% (based on 243 observations) to 18.3% (based on 146 observations) (see Appendix A.1.5). After reducing the sample size to approximately 95 (0.4 birch trees per plot), BP lost significance (see Figure A7). The BP of spruce regeneration decreased from 2015 (12,280 observations) to 2018 (10,774 observations) by −0.5 PP (see Appendix A.1.1). We ware able to eliminate observations to a sample size of 5020 without losing significance (see Figure A8). The change of BP of beech (−6.3 PP between 2015 (4615 seedlings) and 2018 (5250 seedlings)) was measured as significant (Appendix A.1.1). With our sample size reduction, it was only possible to lose significance twice; significance was lost once with a sample size of 895 and once with a sample size of 775 before aborting our computing processes after the error messages reported above (see Figure A8). The occurrence of the two significance losses can be explained by a random combination of the removed sample points and can be interpreted as outliers.
Moreover, for silver fir, we found a significant BP change of +4,6 PP between 2015 (1053 individuals) and 2018 (1062 individuals) (see Appendix A.1.4). The significance remained at a sample size for 679 observations. Thereafter, the lower half of the confidence interval regularly crosses the significance boundary (see Figure A8). Figure A9. Overview of the administrative organization of the study area (Bavarian Forest National Park); RLG = Rachel-Lusen-Region (light blue); FRG = Falkenstein-Rachel-Region (light green); systematic distributed plots of the regeneration inventory. The few plots which fall out of that grid near the NP's borders are not incorporated in our analysis. Figure A10. Spatial distribution of the silver fir (A) and the rowan (B) in our study area.