Does Ungulate Herbivory Translate into Diversity of Woody Plants? A Long-Term Study in a Montane Forest Ecosystem in Austria

Abstract

Different species-specific traits of woody plant species, feeding preferences of herbivores together with resulting effects on plant competition are expected to translate into different plant community structures and expressions of biodiversity. We studied the diversity of woody plant species (trees and shrubs) and structural diversity of forest trees, using a 30-year and an 18-year dataset of ungulate exclosure-control plot pairs in a mixed alpine forest community in Austria. We surveyed the tallest individuals per tree species and plot and analyzed the collective of top-height individuals per plot pair. Incidence data for exclosure and control plots were aggregated. Comparing species diversity and diversity of height classes on the plots throughout time, we calculated diversity profiles based on Hill numbers. Diversity of top height individuals and structural diversity, expressed by height classes, were two diversity aspects that differed between exclosures and control plots. Other diversity estimates of woody plant species showed huge variation without significant differences between plots. Height growth was significantly suppressed by ungulate herbivory. Effects of ungulate herbivores in forest ecosystems are highly complex and context-dependent and thus not reducible to simple top-down forces. Long-term surveys provide data that reflect “ultimate” effects of herbivory interacting with other drivers of community dynamics.

1. Introduction

Wildlife ungulate species are essential drivers of the dynamics of many terrestrial ecosystems, with impacts on biodiversity at several scales. Herbivory affects species-specific mortality patterns of plants and the related interspecific competition [1,2]. Ungulate browsing induces defense responses in plants and modifies their chemical status or their morphology [3]. These defense mechanisms reflect long-term effects of herbivory on plant evolution [4]. Furthermore, ungulate herbivory might also significantly change habitat conditions for other animal taxa, such as birds or arthropods, translating into modified species composition, abundance, biomass, and diversity [5,6,7,8].

For woody plant species, feeding preferences of herbivores are linked to several crucial traits of the plants such as their nutritional value, their taste and palatability, or toxin content [9,10]. The latter require detoxification mechanisms [11] and cause detoxification costs in herbivores [12,13,14]. The combination of these traits might explain general patterns of preferred vs. less preferred woody plant species for browsing. Species-specific traits of plants also include their tolerance towards browsing, ranging from highly sensitive species, which disappear in the presence of specific herbivores, to highly tolerant plant species, which might survive and grow even when being subjected to high browsing pressure and which, in turn, might outcompete more sensitive species [15,16,17]. Selective herbivory [18] thus can induce a loss of preferred plant species and might favor non-preferred and/or tolerant tree species, both potentially translating into different biodiversity measures.

Early life stages of trees are more prone to die from herbivory [19]. Older stages mainly provide nutrient-poor, but carbon-rich tissues to ungulate herbivores [20,21] at those height strata that can be accessed by ungulates. In parallel, particularly ungulates show preferred heights for browsing within their reach [22], which might be related to trade-offs in terms of food choice, biokinetic factors, and predator vigilance [23,24], but also plant morphological traits [25,26]. Aside from this, the accessibility of plant biomass for ungulates is shaped by snow conditions in winter: Early life stages of plants might be covered by snow and thus become inaccessible, while higher plants might be accessed at higher vertical levels due to the prevailing snow heights [27]. Thus, ungulate browsing does not only impact the occurrence of woody plant species and their inter-specific competition, but also the vertical structure and layering of woody plant communities.

Considering the longevity and reachable height of woody plants and trees in particular, ungulate herbivory only covers a comparable small timeframe, which nonetheless might have significant impact on later developmental stages. However, other drivers apart from herbivory like microsite conditions, pathogens, or intra- and interspecific plant competition shape community structures of woody plants at early life stages [28,29]. Consequently, ungulate herbivory not necessarily explains a high proportion of variation in tree growth performance [30] but constitutes one of various interacting impact factors [31]. In turn, dynamics of tree regeneration in terms of species occurrence or height growth patterns distinctly shift over time. As a consequence, results from surveys with short and medium duration differ from long-term results, when individuals have reached later developmental stages. At these later stages, the terminal heights are no longer accessible for browsing by ungulates as is the case with thickets up to mature stands. Initial herbivore effects at early developmental stages might be accompanied by a multitude of co-occurring or subsequent processes, and ungulate herbivory might not be a significant driver of dynamics at later stages, but only form a part of compensatory mortality [30]. As a consequence, short-time surveys on the effects of herbivory, quantifying reductions in stem numbers or net primary production, might overestimate these effects on the long run and at the community level. This is particularly relevant for tree species composition and a well-known drawback in studies on tree regeneration in forest gaps. It has been shown that species composition at early stages (expressed by seedling densities) might not be indicative of the future species composition of forest stands [32,33,34]. Thus, tracing community patterns and the contribution of ungulate herbivory needs long-term surveys under different site conditions [30,35].

Responses of plant communities to mammalian herbivory are typically investigated with exclosure experiments [35,36,37,38], where plant densities or species occurrences within a fence (mammals/ungulates excluded) are contrasted against an open control plot that mammals have access to. Comparing exclosures and coupled open access plots, the consumption of organic matter by herbivores might be expressed as differences in net primary production [39], in plant densities [40,41,42], or, at the community level, in structural diversity or species diversity [17,43,44,45,46].

To compare different aspects and indicators of biodiversity at the community level (species and structural diversity) and potential impacts of ungulate herbivory, we used a unique, long-term data set of ungulate exclosures gained in an alpine mixed forest community in Austria. For each tree species, heights of the highest tree individuals (n_max = 6) were recorded in three-year periods together with coverage and heights of shrub species. This allowed us to address the following questions: (1) Is the diversity of tree species, represented by the highest individuals per species, shaped by ungulate herbivory? (2) How does herbivory affect the diversity of tree species within the class of dominant, top-height individuals per plot pair? (3) Is the structural diversity, expressed by the number of observed tree height classes, impacted by ungulate herbivory? (4) Does shrub diversity in exclosures deviate from the one on open-access control plots? (5) Do differences in diversity change over time?

2. Materials and Methods

2.1. Study Area

Regeneration dynamics were monitored on the mountain massif “Höllengebirge” (highest peak at 1862 m), which comprises a forest area of 90 km². It is part of the northern Limestone Alps in Upper Austria, Central Europe (47°49′ N, 13°30′ E). The forest encloses a central plateau laying above the forest line. The dominating geological substrates are limestone and dolomite, with typical soils being classified as Rendzinas and relictic loams. The climate of the region is submaritime, with long winter periods and cool, wet, and short summer periods. Mean annual temperature of the area ranges between 6 and 8 °C (depending on altitude) and mean annual precipitation reaches approx. 1800 mm. Annual precipitation is characterized by a bimodal temporal pattern, with one maximum in July and another one between November and January. The duration of snow cover averages 170 days per year, with large variation between single years. The typical montane mixed forest communities comprise Norway spruce (Picea abies), silver fir (Abies alba), European beech (Fagus sylvatica), European ash (Fraxinus excelsior), and sycamore maple (Acer pseudoplatanus). They are classified as Helleboro nigri-Fagetum and Adenostylo glabrea-Fagetum. Forests are mainly managed via two basic silvicultural systems, one being a permanent forest system with shelterwood and strip-selection cutting and natural forest regeneration, and the second being strip clear-cuts with supplementary afforestation. The investigation is focused on areas with natural forest regeneration.

The estimated mean ungulate density on the forested area in spring (without newborn calves, fawns, and kids) totals about 15–24 individuals per km². It is composed of approximately 4 red deer (Cervus elaphus, counted at winter-feeding stations), 5–10 roe deer (Capreolus capreolus, estimated by local hunters), and 6–10 chamois (Rupicapra rupicapra, estimated by local hunters) per km². The hunting season runs from May to December, yielding hunting bags of 1–2 red deer, 1–4 roe deer, and 1–2 chamois per km² and year on a constant level over the observation period.

2.2. Sampling Method

We initially established 100 pairs of ungulate exclosures (excluding the occurring red deer, roe deer, and chamois) and unfenced control plots at altitudes between 500–1100 m a.s.l. in 1989, and we completed the survey in 2019. Until this point of time, several fence or control plots had been destroyed or disturbed by natural abiotic disturbance events like storm throw, snow loading, or rock falls. Plots that experienced other disturbances than ungulate herbivory were omitted. Thus, a final number of 14 plot pairs could be surveyed over the entire 30-year time span. The minimum distance between plot pairs was 200 m.

Plots pairs were installed on areas where forest regeneration was in an initial phase; regeneration was just starting (maximum tree height of regeneration approx. 0.5 m) or was soon expected to occur. In choosing sites, it was also considered that light conditions would enable further development of seedlings (maximum canopy density by seed trees about 80%) and that a subsequent dense canopy closure was unlikely. The 14 plot pairs that were investigated for 30 years were regenerated via natural regeneration (without planting). At the beginning of the investigation, the canopy of these plots was less than 50% in 6 cases, and between 50 and 80% in 8 cases. After 30 years, the canopy was reduced in most cases; 12 plot pairs had less than 50% or no canopy anymore, only 2 pairs had more than 50% canopy.

The mature forest stands, in which the plot pairs were situated, were subject to regular forest management, i.e., cutting of mature trees, and successively reducing canopy. On the plots themselves, no forest interventions were carried out, i.e., no thinning of forest regeneration (no cutting of young trees). If the regeneration of a plot pair was influenced unequally between a fenced and paired unfenced plot by forest operations in the mature stand, this plot pair was omitted from the investigation.

Smaller mammals like hares (Lepus europaeus) or herbivore voles and mice were able to enter the exclosures. The type of treatment (fencing vs. control) was randomly attributed with a coin toss. Exclosures were constructed with a 2 m tall, galvanized fence outlining an area of 6 × 6 m. Distances between exclosures and corresponding control plots ranged from 5 m up to max. 20 m, with the largest possible congruence of site and stand conditions between exclosures and control plots. The entire sampling procedure was optimized to capture effects of ungulate herbivory and to keep other driving factors as constant as possible. Fixing of distances between plot pairs would have meant a severe constraint in the search for plots with the highest possible similarity between exclosure and control plot pairs. The distance range of 5 to 20 m between plot pairs reduces mutual impacts of paired plots on the one hand, but also offers a larger search buffer to find comparable conditions for forest regeneration. This is particularly necessary in steep mountainous regions.

In this alpine environment with very high spatial heterogeneity, larger plots would have inevitably implicated more divergence between fence and control plots even at the beginning of the survey, which would have biased our records. Within the 6 × 6 m area of each fence and control plot, we recorded plant data within the central 5 × 5 m square (i.e., survey area without “edge effects”), which was permanently marked with metal rods. Data were first recorded in the year of establishment of exclosure/control pairs and afterwards in a three-year cycle, except for the last phase of the survey, which covered a 6-year time span (summing up nine survey cycles for the longest lasting surveys); in 2016, we did not record the regeneration on the plots.

Considering the gap-maker–gap-filler approach, we dynamically recorded the height of the largest trees per tree species in paired plots with and without exclosures, as this allowed us to assess herbivore impacts on those individuals, which have a high probability of attaining reproductive stages. We thus recorded the six highest individuals per tree species and plot or, in case of lower stem numbers, all occurring individuals per tree species and sampling plot (25 m², each fenced and unfenced) [see 30]. We used the following height classes: ≤10 cm, 11–25 cm, 26–40 cm, 41–70 cm, 71–100 cm, 101–130 cm, 131–160 cm, 161–200 cm, 201–250 cm, 251–300 cm,… (continuing 50 cm classes). These height classes cover an increasing height range with increasing height, starting with 15 cm intervals for seedlings and ending up with 50 cm intervals at later stages. These unequal class ranges were used for two reasons: (1) They account for typical height growth patterns in tree species with steep slopes at early life stages and flattening curves after decades [47]; (2) consequently, small differences in height or height growth might be more meaningful in terms of competitive power and survival in smaller trees compared to later stages [29,48,49]. For shrubs, we visually assessed the coverage (%) of shrub species per survey plot, as individual entities and related numbers might not be easily identified in the course of a field survey.

No animal experiments were carried out. The ungulates were neither caught nor influenced by medication. They were only denied access to small patches within forests by means of a fence. The fencing was carried out in accordance with all relevant guidelines and regulations. The entire sampling design was optimized to capture effects of ungulate herbivory and to keep other driving factors as constant as possible.

2.3. Data Preparation

As the height class distribution of the sampled trees might distinctly deviate even within the collective of the highest six individuals per tree species, we additionally defined a group of “top height individuals”, which reach or exceed the 66%-height-level of the highest tree individual per exclosure/control pair. The threshold of 66% was set arbitrarily, representing the upper third of the height range. This subsample of tree individuals represents the collective of dominant stems (hereinafter termed “top height individuals”) with presumed highest competitive power. Due to a delay in height growth between some plot pairs, tree individuals of one plot did not reach the 66% threshold of the maximum value of the paired plot in several cases. These plot pairs had to be omitted from the analysis of top height individuals.

Addressing structural diversity of tree species, we looked at the occurrence of different height classes in the exclosures and on the corresponding open control plots. As height classes cover smaller height ranges at the seedling stage and change to larger range values for higher tree individuals, this height class diversity does not reflect absolute differences in height growth between the exclosure/control pairs, but depicts potential differences in terms of observed height class numbers.

To additionally cover the issue of potentially impacted height growth patterns of woody plants due to ungulate herbivory, we also analyzed height trajectories of the highest six individuals per tree species. As tree heights had been recorded in terms of height classes, we used the mid-point height per height class interval as representative per individual tree (e.g., 85 cm for the 71–100 cm height class).

For shrubs, a shrub coverage above zero was interpreted as presence of a given shrub species.

2.4. Diversity Indices and Statistical Analyses

Woody plants (trees and shrubs) of the same species are often spatially aggregated, and individuals should not be interpreted or modeled as independent sampling units [50]. Accounting for this potential impact on woody plant species distribution and assemblages, we converted the abundance data of tree and shrub species and height classes of trees in each study plot into incidence-based (detection versus non-detection) data. Study plots were afterwards treated as sampling units, and incidence data for exclosure and control plots were aggregated and analyzed based on a species-by-sampling-unit incidence matrix [50,51]. Analogous to an abundance-based model, the proportion of study plots in which a species can be detected is represented by the incidence probability [50].

To compare woody plant species diversity and the diversity of height classes of trees of exclosure and control plots and throughout time, we calculated diversity estimates based on the effective number of species. The first three Hill numbers [52] of order q thereby refer to species richness (q = 0, hereinafter referred to as “species richness”), the exponential of Shannon’s entropy index (q = 1, hereinafter referred to as “Shannon diversity”), and the inverse Simpson´s concentration index (q = 2, hereinafter referred to as “Simpson diversity”) [51].

Differences in diversity estimates in terms of the three Hill numbers can be briefly summarized, based on detailed explanations described in several publications [50,51,52]: Diversity estimates based on a Hill number of order q = 0 are highly sensitive regarding infrequent species, because all species (frequent and infrequent) are counted equally, without accounting for their detection probability. For estimation of diversity based on a Hill number of order q = 1, species are weighted by their probability of occurrence within the sample. This probability is affected by the area sampled and the rarity of each species within the study area. Thus, diversity estimates based on a Hill number of order q = 1 can be interpreted as effective number of frequent species in the assemblage. Diversity estimates based on a Hill number of order q = 2 are highly sensitive to highly frequent species, because species are proportionally weighted by its squared detection probability. Therefore, diversity estimates of order q = 2 can be interpreted as the effective number of highly frequent species in the assemblage.

We then followed a procedure described in detail by Colwell et al. (2012) [53] and Chao et al. (2020) [50]: We calculated sample size-based rarefaction and extrapolation curves, as detection or occurrence of species is highly dependent on the number of sampling units (study plots). The sample size was extrapolated to double the size of the surveyed exclosure or control plots [51]. It is assumed that the observed diversity is equivalent to the “true” diversity if the sample size-based rarefaction and extrapolation curve stabilizes and approaches an asymptote, thereby enabling comparisons of diversities across multiple assemblages (exclosure versus control plots) [50]. Sufficiency of data (sample completeness) was determined by visually examining empirical diversity profiles of sample size-based rarefaction and extrapolation curves for orders q. Sample size-based rarefaction and extrapolation curves usually level off for Hill numbers of order q = 1 (Shannon diversity) and q = 2 (Simpson diversity), thereby indicating that the sample size is sufficient for interpreting frequent and highly frequent species in the assemblages, respectively. In contrast, when sample size-based rarefaction and extrapolation curves do not stabilize (typically for Hill numbers of order q = 0, species richness), the asymptotic estimator represents only a lower bound, because some rare species might not occur within the set of study plots. In such cases, a non-asymptotic standardization approach was used [50]. We then compared diversity estimates of exclosure and control plots by calculating diversity estimates based on a standardized sample coverage (maximum coverage C_max) [50]. Overall, 95% confidence intervals were obtained from a bootstrap method [50,51,54]. We performed preliminary simulations with varying numbers of replications for the bootstrap method and our simulations suggested that stable variance estimates and confidence intervals were obtained for our dataset when using 400 replications. Non-overlap of confidence intervals indicated that diversity measures of exclosure and control plots differ significantly (p < 0.05) [53]. Accounting for the restricted number of plots pairs that could be surveyed over the entire 30-year period (n = 14), we also ran calculations for the plots that could be monitored for 18 years (n = 43). This was particularly meaningful for those analyses, where we defined a specific subsample (e.g., the top-height individuals). In cases where the two sample sizes did not yield different results for study year 18, we report the outcomes for the 30-year period and provide results of the 18-year period in the Supplementary material. In case that the larger sample size at year 18 yielded different results and smaller sample sizes did not produce stables variance estimates and confidence intervals, we only present the calculations for the larger sample (18 years).

We tested for differences in tree height trajectories within the 30-year study period and between exclosure and control plots by using a Kruskal-Wallis rank sum test followed by Wilcoxon rank sum tests adjusted for multiple comparisons using Bonferroni correction (significance level of adjusted Wilcoxon rank sum tests p < 0.05). All analyses were performed using the R packages ‘iNEXT’ [55] and ‘iNEXT.4steps’ [50]. Results were visualized using the R packages ‘ggplot2′ [56] and ‘ggpubr’ [57].

3. Results

3.1. Tree Species Diversity

Within the 14 exclosure-control pairs, tree individuals completely disappeared on one control plot 6 years after the establishment of the plots. This means that no trees survived until the last survey due to ungulate impacts contrary to the fenced plot, while all other conditions remained comparable. The zero-value of tree species incidence on the open plot did not allow for further statistical calculations. Nevertheless, it was not a bias that would have led to an a priori exclusion of the respective fence-control pair, but trees did not survive on this plot.

A total of 14 different tree species were detected on the remaining 26 study plots (13 exclosure and 13 control plots) within the 30-year study period (see

Table 1. Descriptive statistics (min, max, median) of species richness (number of tree species detected) within the 13 control and exclosure plot pairs, depicted for three points of time within the 30-year study period.

Year	Control Plots				Exclosure Plots
	Min	Max	Median	n	Min	Max	Median	n
9	1	7	5	13	3	9	6	13
18	2	7	4	13	3	8	6	13
30	3	7	4	13	3	7	5	13

). On the 13 control plots, the observed number of tree species varied between 1 and 7 species over the whole study period. In contrast, we recorded 3 to 9 species on the exclosure plots in the same time span. A species list is available in the Supplementary material (Supplementary Material Table S1).

Sample size-based rarefaction and extrapolation sampling curves for Shannon diversity (q = 1) and Simpson diversity (q = 2) levelled off and approached an asymptote, thereby indicating that the number of sampling plots was sufficient to represent frequent and highly frequent tree species in the assemblages (Figure 1). For species richness (q = 0), the increasing sample size-based rarefaction and extrapolation curves for exclosure plots revealed that data were insufficient to infer “true” species richness, because some rare tree species might not have been detected within the set of study plots but could have occurred in the surrounding study region (Figure 1, left panel). Thus, the asymptotic estimator represented only a lower bound and rare tree species might have been underrepresented. Testing for differences in species richness between exclosure and control plots was therefore made at a standardized coverage value of C_max = 96.3% (

Table 2. Asymptotic and non-asymptotic diversity estimates for incidence-based tree species data collected on control and exclosure plots in Austria 9, 18, and 30 years after establishment of study plots. Non-asymptotic rarefaction and extrapolation analyses are based on point diversities for a specified sample coverage of C_max = 96.3%. Numerical values for Hill numbers of order q (q = 0, species richness; q = 1, Shannon diversity; q = 2, Simpson diversity).

	Control Plots			Exclosure Plots
Year	q = 0	q = 1	q = 2	q = 0	q = 1	q = 2
Asymptotic diversity estimate (mean ± SE)
9	10.92 ± 1.21	8.85 ± 0.52	8.14 ± 0.52	16.15 ± 4.43	9.94 ± 0.85	8.51 ± 0.53
18	9.46 ± 1.41	8.01 ± 0.56	7.36 ± 0.49	22.23 ± 6.29	10.53 ± 1.06	8.66 ± 0.55
30	9.46 ± 1.42	7.97 ± 0.51	7.36 ± 0.43	15.15 ± 4.35	9.02 ± 0.86	7.76 ± 0.49
Non-asymptotic point diversity estimate [95% confidence interval: lower; upper]
9	9.31 [7.5; 11.1]	8.09 [7.1; 9.1]	7.49 [6.6; 8.3]	12.13 [6.5; 17.7]	9.24 [7.5; 10.9]	8.22 [7.2; 9.3]
18	8.29 [6.3; 10.3]	7.29 [6.2; 8.4]	6.72 [5.8; 7.6]	16.80 [7.5; 26.2]	10.00 [7.7; 12.3]	8.51 [7.3; 9.7]
30	8.19 [6.3; 10.0]	7.28 [6.3; 8.3]	6.77 [6.0; 7.6]	11.35 [5.5; 17.2]	8.44 [6.8; 10.1]	7.54 [6.5; 8.6]

, Figure S1). Overall, 95% confidence intervals of species richness (Figure S1), as well as Shannon diversity and Simpson diversity (Figure 1) of tree species on exclosure plots and control plots, largely overlapped at the three points of time, explicitly addressed in our analyses (Table 2). Results for tree species diversity within the larger sample size of the 18-year period are shown in the Supplementary material (Table S3, Figures S2 and S3).

3.2. Tree Species Diversity of Top Height Individuals

Due to a delay in height growth, tree individuals on several plots did not reach the 66% height-threshold of the paired plot. Thus, these plot pairs could not be included in the calculations and the sample size was therefore reduced. This fact led to a sample size for the 30-year period (n_{plot pairs} = 7), which was too small to produce stable confidence intervals. Thus, we used the 18-year samples (n = 22) for the analysis of the tree species diversity of top height individuals.

A total of 14 different tree species were detected on the 22 study plot pairs for the collective of top height individuals 18 years after establishment of study plots. On the control plots, the observed number of tree species varied between 1 and 4 (median = 4) species, compared to 1 up to 7 (median = 3) observed tree species on the exclosure plots.

The number of sampling plots was sufficient to represent frequent (q = 1) and highly frequent (q = 2) species in the assemblages, because sampling curves approached an asymptote (Figure 2).

In the survey that was conducted 18 years after establishment of study plots, Shannon (q = 1) and Simpson (q = 2) diversity estimates showed that tree species diversity of top height individuals was significantly higher in exclosure plots in comparison to control plots. Based on the asymptotic results, the difference in the observed tree species diversity between exclosure and control plots was 6.40 with respect to frequent tree species (q = 1) and 5.43 with respect to highly frequent tree species (q = 2) 18 years after establishment of study plots (

Table 3. Asymptotic and non-asymptotic diversity estimates for incidence-based tree species data of top height individuals collected on control and exclosure plots in Austria 18 years after establishment of study plots. Non-asymptotic rarefaction and extrapolation analyses are based on point diversities for a specified sample coverage of Cmax = 98.8%. Numerical values for Hill numbers of order q (q = 0, species richness; q = 1, Shannon diversity; q = 2, Simpson diversity).

	Control Plots			Exclosure Plots
Year	q = 0	q = 1	q = 2	q = 0	q = 1	q = 2
Asymptotic diversity estimate (mean ± SE)
18	6.24 ± 1.85	4.22 ± 0.64	2.98 ± 0.48	16.55 ± 5.41	10.62 ± 1.08	8.41 ± 0.82
Non-asymptotic point diversity estimate [95% confidence interval: lower; upper]
18	6.14 [2.6; 9.6]	3.94 [2.8; 5.1]	2.86 [2.0; 3.7]	15.99 [3.4; 28.6]	10.22 [7.9; 12.5]	7.11 [6.6; 9.6]

, Figure 2).

The increasing sample size-based rarefaction and extrapolation curves for study plots surveyed 18 years after plots were established revealed that data were insufficient to infer “true” tree species richness (q = 0), because some rare tree species might have not been detected within the set of study plots, but could have occurred in the study region (Figure 2). Thus, we used a standardized coverage value of C_max = 98.8% to test for differences in tree species richness between exclosure and control plots. Thereby, 95% confidence intervals of tree species richness estimates (q = 0) in exclosure plots and control plots overlapped in study plots surveyed 18 years after implementation (Table 3, Figure 3).

3.3. Height Class Diversity of Tree Species

At the three notably timepoints in our 30-year-long study, particularly addressed in the given analyses (i.e., year 9, year 18, and year 30), 35 different tree height classes were recorded on the study plots (see

Table 4. Descriptive statistics (min, max, median) of height class richness (number of tree height classes detected) on the 13 control and 13 exclosure plots, depicted for three points of time within the 30-year study period.

Year	Control Plots				Exclosure Plots
	Min	Max	Median	n	Min	Max	Median	n
9	2	7	4	13	3	7	5	13
18	2	8	4	13	4	10	7	13
30	4	10	6	13	4	13	9	13

). Numbers of height classes ranged between 2 and 10 on the exclosure plots, whereas 3 to 13 height classes were observed on the control plots.

Sample size-based rarefaction and extrapolation sampling curves with 95% confidence intervals for Shannon diversity (q = 1) and Simpson diversity (q = 2) levelled off and approached an asymptote, thereby indicating that the number of sampling plots was sufficient to represent frequent and highly frequent height classes of trees in the assemblages (Figure 4). In surveys, that were conducted 9, 18, and 30 years after establishment of study plots, Shannon (q = 1) and Simpson (q = 2) diversity estimates of tree height classes showed that tree height class diversity was significantly higher in exclosure plots in comparison to control plots. Based on the asymptotic results, the difference in observed tree height classes between exclosure and control plots was 2.20, 3.98, and 10.63 with respect to frequent height classes (q = 1) 9, 18, and 30 years after establishment of study plots, respectively (

Table 5. Asymptotic and non-asymptotic diversity estimates for incidence-based tree height class data collected on control and exclosure plots in Austria 9, 18, and 30 years after establishment of study plots. Non-asymptotic rarefaction and extrapolation analyses are based on point diversities for a specified sample coverage of C_max = 96.4%. Numerical values for Hill numbers of order q (q = 0, height class richness; q = 1, Shannon diversity; q = 2, Simpson diversity).

	Control Plots			Exclosure Plots
Year	q = 0	q = 1	q = 2	q = 0	q = 1	q = 2
Asymptotic diversity estimate (mean ± SE)
9	8.92 ± 1.46	6.34 ± 0.48	5.48 ± 0.39	11.85 ± 2.31	8.54 ± 0.66	7.65 ± 0.47
18	14.77 ± 3.44	10.49 ± 0.99	8.87 ± 0.83	20.15 ± 4.76	14.47 ± 1.07	13.20 ± 0.76
30	30.38 ± 10.04	20.51 ± 2.14	15.23 ± 1.60	44.54 ± 10.00	31.14 ± 2.54	25.17 ± 1.99
Non-asymptotic point diversity estimate [95% confidence interval: lower; upper]
9	7.65 [5.8; 9.5]	5.96 [5.1; 6.8]	5.33 [4.6; 6.0]	9.40 [6.4; 12.4]	7.92 [6.6; 9.2]	7.27 [6.3; 8.2]
18	12.67 [8.0; 17.4]	9.71 [7.9; 11.5]	8.37 [7.0; 9.7]	15.27 [9.7; 20.8]	13.14 [11.0; 15.3]	12.15 [10.6; 13.7]
30	27.64 [14.2; 41.0]	19.03 [15.1; 23.0]	14.33 [11.8; 16.8]	39.55 [21.0; 58.1]	29.00 [23.7; 34.3]	23.44 [20.0; 26.8]

, Figure 4). The difference with respect to highly frequent tree height classes (q = 2) between exclosure and control plots increased from 2.17 to 4.33 and 9.94 in the 9, 18, and 30 years after establishment of study plots (Table 5, Figure 4).

The increasing sample size-based rarefaction and extrapolation curves for study plots surveyed 9, 18, and 30 years after plots were implemented revealed that data were insufficient to infer “true” height class richness (q = 0), because some rare height classes might have not been detected within the set of study plots but could have occurred in the study region (Figure 4). Thus, we used a standardized coverage value of C_max = 96.4% to test for differences in height class richness between exclosure and control plots. The 95% confidence intervals of height class richness estimates (q = 0) largely overlapped between exclosure plots and control plots in study plots surveyed 9, 18, and 30 years after implementation (Table 5, Figure S4). Results for height class diversity within the larger sample size of the 18-year period are shown in the Supplementary material (Table S4, Figures S5 and S6).

3.4. Shrub Species Diversity

A total of 19 different shrub species were detected in the study plots within the three consecutive surveys carried out within the 30-year study period after study sites had been implemented (see

Table 6. Descriptive statistics (min, max, median) of species richness (number of shrub species detected) within the control and exclosure plots throughout the 30-year study period.

Year	Control Plots				Exclosure Plots
	Min	Max	Median	n	Min	Max	Median	n
9	1	4	2	13	1	4	3	13
18	1	8	2	13	1	4	2	13
30	1	8	2	11	1	7	2	11

). On two plot pairs, where trees reached maximum heights of up to 18 m, shrubs disappeared both on the fenced and the control plot, reducing the sample size from 13 to 11 plot pairs for the last survey period (Table 6). In the control plots, the observed number of shrub species varied between 1 and 8 species throughout the whole study period. In contrast, the number of observed shrub species in exclosure plots varied between 1 and 7 species in the same period. A species list is available in the Supplementary material (Table S2).

Increasing sample size-based rarefaction and extrapolation curves for exclosure study plots revealed that data were insufficient to infer “true” species richness (q = 0), because some rare shrub species might have not been detected within the set of study plots but could have occurred in the study region (Figure 5). Testing for differences in species richness between exclosure and control plots was therefore made at a standardized coverage value of C_max = 90.1%. In contrast, sample size-based rarefaction and extrapolation sampling curves for Shannon diversity (q = 1) and Simpson diversity (q = 2) levelled off and approached an asymptote, thereby indicating that the number of sampling plots was sufficient to represent frequent and highly frequent shrub species in the assemblages (Figure 5). Overall, 95% confidence intervals of species richness (Figure S7), Shannon diversity, and Simpson diversity (Figure 5) of shrub species in exclosure plots and control plots largely overlapped (

Table 7. Asymptotic and non-asymptotic diversity estimates for incidence-based shrub species data collected within control and exclosure plots in Austria 9, 18, and 30 years after establishment of study plots. Non-asymptotic rarefaction and extrapolation analyses are based on point diversities for a specified sample coverage of C_max = 90.1%. Numerical values for Hill numbers of order q (q = 0, height class richness; q = 1, Shannon diversity; q = 2, Simpson diversity).

	Control Plots			Exclosure Plots
Year	q = 0	q = 1	q = 2	q = 0	q = 1	q = 2
Asymptotic diversity estimate (mean ± SE)
9	16.77 ± 7.81	11.14 ± 2.31	8.26 ± 1.54	14.15 ± 4.68	10.03 ± 1.58	8.89 ± 1.21
18	17.77 ± 7.51	12.61 ± 2.40	9.86 ± 1.89	13.69 ± 5.52	9.72 ± 1.79	7.76 ± 1.35
30	25.70 ± 11.90	18.65 ± 3.78	13.47 ± 2.65	14.84 ± 6.21	12.20 ± 2.02	9.80 ± 1.63
Non-asymptotic point diversity estimate [95% confidence interval: lower; upper]
9	13.24 [2.8; 23.6]	9.41 [5.0; 13.8]	7.30 [4.7; 9.9]	9.27 [4.2; 14.4]	7.94 [5.3; 10.6]	7.17 [5.2; 9.1]
18	14.00 [4.1; 23.9]	10.56 [5.8; 15.3]	8.50 [5.5; 11.5]	10.83 [3.0; 18.7]	8.12 [4.8, 11.4]	6.73 [4.7; 8.8]
30	21.32 [4.5; 38.1]	15.87 [8.7; 23.1]	11.69 [7.5; 15.9]	12.95 [4.3, 21.7]	10.07 [6.4; 13.8]	8.18 [5.7; 10.6]

, Figure 5). Results for shrub species diversity within the larger sample size of the 18-year period are shown in the Supplementary material (Table S5, Figures S8 and S9).

3.5. Tree Height

Heights of trees varied throughout time and between exclosure and control plots (Kruskal-Wallis rank sum test, Chi-squared = 452.42, df = 5, p < 0.001, see Figure 6). Within control and exclosure plots, tree height significantly increased throughout the 30-year study period (Figure 6, Supplementary material Table S6). In addition, trees were significantly higher on the exclosure plots in comparison to control plots 9, 18, and 30 years after study plots were implemented (Figure 6,

Table 8. Descriptive statistics (min, max, median) of height (cm) of the tallest individuals (n_max = 6) per tree species detected on the 13 control and 13 exclosure plots throughout the 30-year study period and result of pairwise comparisons using Wilcoxon rank sum test with Bonferroni correction. Significance level of adjusted Wilcoxon rank sum tests p < 0.05.

Year	Control Plots				Exclosure Plots
	Min	Max	Median	n	Min	Max	Median	n	p-Value
9	5	225	18	259	5	325	33	292	<0.001
18	5	475	33	225	5	725	145	265	<0.001
30	5	1575	55	168	5	1825	180	210	<0.001

).

4. Discussion

Herbivory is seen as a key process in forest ecosystems that influences species occurrence and ecosystem functions [58,59,60,61]. Feeding of large herbivores can have effects on ecosystems´ diversity at different levels: Plant species that are highly attractive for (specific guilds of) herbivores in terms of nutritive value, palatability, detoxification costs, or availability might be intensely browsed by plant feeders [9,13,62,63] and thus be less competitive or even inferior compared to plants with other traits [64]. However, the ultimate effects of these plant traits also strongly depend on the physiology and specialization of the occurring or dominating herbivore species [65] and their resulting preference for certain plant species [18,64,66]. Availability of alternative feeding resources, their amount and distribution [67,68], together with herbivore densities can be decisive for the survival and occurrence of such species. On the other hand, plant species develop a series of specific strategies to escape from ungulate herbivory, e.g., by optimizing resistance and tolerance [10]. In trees, species-specific and context-dependent growth patterns [10] also determine the time span, where palatable plant parts are accessible for mammalian herbivores. Further, highly preferred plant species might detract herbivores from other plant species, which thereby escape feeding pressure. For ungulate herbivores, rare plant species can be highly attractive for browsing [69]. Ungulate herbivory might not only impact the occurrence and abundance of plant species by feeding on available and preferred plants, but ungulates might also maintain a higher species richness of specific plant strata [46]. This multitude of functional responses might explain the variability of outcomes from studies that address herbivores’ effects on plant species diversity in forest ecosystems, ranging from decreased to enhanced diversity [e.g., 35]. Our results are in line with this overall picture of high variability in terms of tree species diversity, particularly expressed by huge confidence intervals of our richness profiles. Moreover, the largely overlapping confidence intervals of all three applied diversity measures (species richness, Shannon diversity, and Simpson diversity) at any point of time clearly reflect the versatile responses of tree species in terms of diversity without a clear differentiation between the collective of open plots and exclosures. For Shannon diversity and Simpson diversity, the side-by-side position of diversity profiles of open plots and exclosures highlights the lack of a clear and consistent impact of ungulate herbivory on tree species diversity. The overlapping confidence intervals indicate that impacts of herbivory on tree species diversity might both include suppressed diversity due to ungulate herbivory on some subgroups of individual plot pairs, but also enhanced diversity on other plots.

Within a patch of forest regeneration, the highest, dominating tree individuals might be particularly exposed to browsing at early life stages due to preferred foraging heights of ungulates [22], or a lack of protecting snow cover in winter. It might thus be expected that these dominant tree individuals respond distinctly to an exclosure of ungulate herbivores in terms of changed species diversity. We defined the group of “top-height individuals” as the collective of stems, that reached or exceeded the 66%-height-level of the highest tree individual per exclosure/control pair. This subsample of tree individuals represents the collective of dominant stems with presumed highest competitive power and the highest chance of ingrowth into later developmental stages. For the 18-year period, diversity of frequent and highly frequent tree species within this collective of top height individuals was significantly higher in the exclosures. However, the diversity estimates for species richness (based on the standardized sample coverage) did not differ in terms of 95% confidence intervals. This clearly implicates that ungulate herbivory did translate into a distinct modification of species diversity of top height trees for the collective of fenced and control plots in our study region.

In northern temperate forest ecosystems, snow cover and the resulting availability of plant parts determine the net biomass that is consumed by herbivores [27,70]. Compared to other plant strata and taxa, woody plants (and particularly deciduous species) on the one hand show distinct seasonal changes in the amount of biomass being usable for large herbivores. On the other hand, they constitute the most important feeding resource during periods with snow cover [71,72]. However, co-occurring woody plant species are potential competitors and intensive herbivore feeding on one species or specific individuals might shape inter- and intraspecific competition for the others [28]. For ungulate browsers, a “typical” order of decreasing preference of plant strata has been provided by several authors [73,74,75], starting with the highest preference of grasses, followed by forbs and shrubs, and ending with trees. For sure, these preferences depend on the specific, prevailing ecosystem and the basic foraging types of occurring ungulate species (see [76]), but they also show that more distinct effects of herbivory on the diversity of shrub species might be expected compared to tree species. Nonetheless, we did not find any indication for significant modification of shrub species diversity, expressed by the three diversity indices. Contrary to the diversity profiles of trees, the diversity curves of open control plots were located above the curve of exclosure plots within the shrub species. This observation is to some degree in line with findings of Pekin et al. (2014) [77] in seasonally arid conifer forests, where ungulate browsing could suppress dominant shrub species in favor of shrub species diversity. Contrarily, Reimoser et al. (2022) [35] observed a slightly negative impact of ungulate browsing on shrub species diversity for a 20-year observation period in the Danube flood plains, Austria. Apart from the feeding preference of ungulate browsers and their potential impact on woody plant species, shrubs might also be outcompeted by faster growing trees. However, the close location of diversity curves for shrub species in our study together with largely overlapping confidence intervals did not allow for any inference on either enhanced or depressed diversity of shrub species by ungulate herbivory at the community level.

Particularly ungulates show preferred plant heights for browsing, which allow for better vigilance, meet biokinetic demands [23,24], and respond to plant morphological traits [25,26]. On the other hand, plants try to vertically escape ungulate herbivory by height increment, by changing leaf traits (photosynthesis, morphology, and chemistry) at specific height strata [10,78], and plant architecture [79]. Together with the preference or avoidance of specific plant species, the duration and intensity of browsing and the seasonally changing food availability, these height preferences might yield modifications of structural features of plant communities. However, the multitude of co-occurring processes that drive plant survival and shape community structure [e.g., 48] explains why only a comparatively small proportion of height growth variance might be traced back to ungulate herbivory in documented cases [30]. Considering these potential pathways, we analyzed impacts of ungulate herbivory on height growth performance and on the related structural diversity, represented by the number of occurring height classes. This aspect of vertical structuring yielded diversity indices, where we found distinct differences between the open control plots and the exclosure plots with higher structural diversity on exclosure plots in terms of Shannon diversity and Simpson diversity. This indicates that frequent and highly frequent tree height classes are more sensitive than infrequent height classes with respect to herbivory. For the richness of height classes, large and overlapping confidence intervals again did not support any differences. Subsuming, we found indications for lower structural diversity due to ungulate browsing in our study context. Thereby, it should be highlighted, that we did not use constant height class ranges in our sampling design for all potentially occurring height classes, but smaller ranges for initial life stages of forest trees and higher ranges for later stages. Thus, height class diversity can be interpreted as a proxy for vertical structural diversity, but not for conclusions on absolute height differences or an absolute degree of layering.

To capture potential effects of ungulate browsing on height dynamics of forest trees, we additionally analyzed the heights of the (six) highest sampled tree individuals per tree species and per plot. While browsing of terminal shoots in a single year does not necessarily induce height growth suppression in individual trees, but might even promote height growth, longer lasting herbivory frequently results in a height growth delay of browsed trees (e.g., [15,80]), lowering their competitive power [38]. In our study, browsing of ungulates did not prevent forest trees from growing with time, but we observed a significantly differing height growth performance at year 9, 18, and 30 of the trial. This clearly reflects a suppression of height growth by ungulates, which primarily occurred between the year 9 and 18 of the survey and could not be compensated afterwards.

In our study, we focused on long-term development of forest regeneration by focusing on tree individuals that have a high probability of attaining reproductive stages and covering a 30-year window of time. Height growth was the only addressed feature that distinctly changed within this timeframe; the other aspects showed particularly high variability between plot pairs, but not over larger time periods.

With our analytic approach, we could account for spatial aggregation of tree species at early life stages and thereby contrast 13 out of 14 exclosure-control plot pairs in terms of diversity for 30 years and 42 out of 43 exclosure-control plots for 18 years. However, we could not involve one exclosure-control plot pair, where no single tree individual remained after 6 years on the open plot. Complete regeneration failure, independently from the type of plot, where it occurs (either exclosure or control plot), cannot be expressed via the chosen diversity measures and is thus not reflected.

Effects of ungulate herbivores in forest ecosystems are highly complex and context-dependent and can thus not be reduced to pure top-down forces [37]. The diversity and inconsistency of study results on ungulate herbivory in forest ecosystems might derive not only from different herbivore–plant community composition and structures [see 35], or from different sampling designs, but also from this complexity, to a large degree. Regeneration dynamics and survivorship of tree species are also the result of neighborhood processes of individuals and species-specific life-history traits together with environmental filters. However, these environmental filters were not the target question we intended to address in this specific study, but their interaction with herbivory could be an interesting focal question of subsequent analyses. Lasting effects of ungulate impacts in early stages of forest regeneration (top-twig browsing, etc.) are only visible decades later because other factors, such as site conditions and forest management, can partly compensate or amplify the former ungulate impacts. Solely modeling forest–ungulate interactions is not sufficient to really understand these complex dynamics and lasting effects. As indicated by our study, it needs robust survey methods feasible in the long run that provide data for developmental phases of forest trees, where considerable parts of tree biomass are not accessible for ungulates anymore and thus reflect “ultimate” effects of herbivory interacting with other drivers of community dynamics.