E ﬀ ects of Twenty Years of Ungulate Browsing on Forest Regeneration at Paneveggio Reserve, Italy

: Forest ecosystems are threatened by di ﬀ erent natural disturbances. Among them, the irruption of large herbivores represents one of the most alarming issues. Several local-scale studies have been carried out to clarify the mechanisms governing ungulate–forest interactions, to understand the e ﬀ ect of wild ungulates overabundance, and to apply conservation plans. However, information at large scales, over long periods of observation and from unmanipulated conditions is still scarce. This study aims to improve our knowledge in this ﬁeld by using repeated inventories to investigate: the types of damage produced by ungulate populations on young trees, the drivers that stimulate browsing activity and its consequences on the speciﬁc composition of seedlings and saplings. To reach these goals, we used data collected during a twenty-year monitoring program (1994–2014) in the forests of Paneveggio-Pale di San Martino Nature Park (Italy). We applied descriptive statistics to summarize the data, GLMs to identify the drivers of browsing activity and Non-Metric Multidimensional Scaling (nMDS) ordinations to investigate the changes in speciﬁc composition of young trees across 20 years. We detected increasing browsing activity from 1994 to 2008 and a decline in 2014. Ungulates browsed preferentially in mature stands, and fed mostly on seedlings and saplings under 150 cm of height. The analysis of the environmental drivers of browsing pressure on the smallest size classes of plants suggests that foraging behavior is inﬂuenced by snowpack conditions, ungulate density and seasonality. Moreover, results underline the fact that ungulates feed mostly on palatable species, especially European rowan, but can also use unpalatable plants as emergency food under high competition levels. nMDS results suggest that rowan seed dispersion might be promoted by deer movements, however, saplings of this species were not able to exceed 30 cm of height because of heavy browsing. This bottleneck e ﬀ ect led to the dominance of unpalatable species, mostly Norway spruce, reducing diversity during forest regeneration. If prolonged, this e ﬀ ect could lead to a reduction of tree species richness, with cascading e ﬀ ects on many parts of the ecosystem, and threatening the resilience of the forest to future disturbances. on plants range browsing This study evidences some of in the on


Introduction
Among herbivores, wild ungulates are among the most troubling species for forest dynamics and forest management, due to their worldwide irruption, the occupation of new ranges by natural  Eiberle and Nigg, 1987 [16].
The percent of browsed trees in size class B was generally higher than size class A. For most species, both size classes showed increasing browsing from 1994 to 2008. This trend stopped in 2014, when the percentage of browsing decreased to levels comparable to those registered in 1994. Some exceptions were European larch, whose increasing trend stopped in 2008, and stone pine and other rare species, whose browsing percentage did not decrease in 2014.

Drivers of Browsing Activity
Among the environmental drivers of browsing probability of trees in size class A, basal area (positive effect) and the percentage of broadleaf trees in the plot (negative effect) showed significant effects (p < 0.01, p = 0.01). For seedlings and saplings in size class B, browsing was negatively influenced by basal area (p < 0.01), tree density (p < 0.01) and aspect (p < 0.01) ( Figure 2).

Experimental Design and Data Collection
Four surveys of forest regeneration phase were conducted in the Park during the summers of 1994, 2003, 2008 and 2014. This sampling frequency was chosen to properly measure the browsing activity and avoid double detections. By spacing sampling seasons by at least five years, we assume that the browsing marks we detected in each survey were caused mostly after the previous survey. The sampling design was based on 101 permanent sampling plots randomly placed along fixed contour lines of 1400, 1600, 1800, 2000 m a.s.l. (see Figure S1B, Supplementary Materials). During the first survey, each plot was permanently marked and site parameters (elevation, aspect, slope) were recorded. Each plot covered an area of 400 m 2 . During each survey, diameter at breast height (DBH) was measured for all trees in each plot with DBH greater than 17.5 cm; common forest structure parameters (tree density, mean diameter, basal area per hectare) were then calculated in the lab. To collect data from young trees, a variable radius sampling was carried out in each plot. Trees smaller than 17.5 cm in DBH were counted within a 12-m radius from the center of each plot, starting from a 0-degree azimuth and moving clockwise until 40 trees were counted. If this number was reached before completing Forests 2020, 11, 612 4 of 12 a full 360 • survey, the total area covered by the measurement was recorded by calculating the area of the circular sector corresponding to the azimuth of the last individual tallied. When less than 40 trees had been counted, the measurement was repeated within a 30 m radius from the plot center. The minimum of 40 juvenile trees sampled in each plot ensured that we could rely on an adequate sample size for individual-level analysis. Plot-level summaries were built based on per-hectare density of juvenile trees, which is invariant relative to plot size. We assumed a homogenous spatial distribution of juvenile trees in each plot.
Trees into the regeneration phase were subdivided in four dimensional classes: class A, from 10 to 30 cm in height; class B, from 30 to 150 cm; class C, taller than 150 cm with DBH smaller than 4 cm; class D, taller than 150 cm with DBH between 4 and 17.5 cm. Seedlings under 10 cm of height were not recorded. Each plant was identified, classified, and all possible different types of ungulate damage (bark stripping, fraying, browsing) were recorded on a binary scale.

Statistical Analyses
To investigate the drivers of browsing activity at an individual level, we performed logistic regression (response variable: Damage/undamaged), using as predictors plant species, elevation, slope, aspect, distance from open areas (determined from existing land cover maps), density of adult trees, basal area and percentage of broadleaf trees (DBH > 17.5 cm) in the plot. Data collected during different sampling seasons were used together to perform this analysis. Aspect values (0 • -360 • ) were transformed using the cosine function to correspond with north (1) and south (−1). The most commonly browsed size classes were A and B, so we built a separate model for each of these size classes, and we excluded C and D classes from this analysis.
To investigate temporal variations in species composition of trees at the regeneration stage, we performed a Non-Metric Multidimensional Scaling (nMDS) using the function metaMDS within the VEGAN package for R. To calculate distances between groups, the Zero Adjusted Bray Curtis index was chosen as the dissimilarity metric. Finally, we used the Analysis of Similarities test (ANOSIM) to check whether we could reject the null hypothesis that the similarity in species composition of young trees between years was greater than or equal to the similarity within the years. All analyses were performed using the R statistical framework (R Core Team, 2020).

Browsing Damage
Among the three types of damage, the most abundant was browsing of apices for all the surveys. Since fraying and bark stripping were rare, we decided to not consider them as a threat for regeneration, and excluded them from all other analyses. Browsing pressure on young trees increased between 1994 and 2008 and decreased in 2014. The most browsed classes were A and B, while C was lightly browsed and D was rarely browsed (Table 1). There was a clear difference among the species in terms of browsing damage. Species that showed the highest browsing percentage were silver fir, European beech and European rowan. For both size classes A and B, browsing pressure on these species were constantly above their respective bearable limit, suggested by Eiberle and Nigg (1987) [16], which indicates the maximum percentage of saplings that can be browsed per year without long-term regeneration failure. The coniferous species stone pine and European larch were lightly browsed and their percentages were around their respective tolerance limits. Norway spruce was also lightly browsed, however the pressure on this species was constantly above its tolerable limit of 12% [16], except for in 2014 ( Figure 1).
The percent of browsed trees in size class B was generally higher than size class A. For most species, both size classes showed increasing browsing from 1994 to 2008. This trend stopped in 2014, when the percentage of browsing decreased to levels comparable to those registered in 1994. Some exceptions were European larch, whose increasing trend stopped in 2008, and stone pine and other rare species, whose browsing percentage did not decrease in 2014.

Drivers of Browsing Activity
Among the environmental drivers of browsing probability of trees in size class A, basal area (positive effect) and the percentage of broadleaf trees in the plot (negative effect) showed significant effects (p < 0.01, p = 0.01). For seedlings and saplings in size class B, browsing was negatively influenced by basal area (p < 0.01), tree density (p < 0.01) and aspect (p < 0.01) ( Figure 2). The percent of browsed trees in size class B was generally higher than size class A. For most species, both size classes showed increasing browsing from 1994 to 2008. This trend stopped in 2014, when the percentage of browsing decreased to levels comparable to those registered in 1994. Some exceptions were European larch, whose increasing trend stopped in 2008, and stone pine and other rare species, whose browsing percentage did not decrease in 2014.

Drivers of Browsing Activity
Among the environmental drivers of browsing probability of trees in size class A, basal area (positive effect) and the percentage of broadleaf trees in the plot (negative effect) showed significant effects (p < 0.01, p = 0.01). For seedlings and saplings in size class B, browsing was negatively influenced by basal area (p < 0.01), tree density (p < 0.01) and aspect (p < 0.01) ( Figure 2). Elevation, slope and distance from open areas did not have any significant effect on the probability of browsing for either size class (see Table S1, Supplementary Materials).
GLMs also confirmed that relative to Picea, which served as reference category, the probability of browsing was significantly higher for Abies, Fagus, Sorbus and other broadleaf trees species, and lower for Pinus trees in class B (Figure 2). Larix and class A Pinus trees were not significantly different. The Elevation, slope and distance from open areas did not have any significant effect on the probability of browsing for either size class (see Table S1, Supplementary Materials).
GLMs also confirmed that relative to Picea, which served as reference category, the probability of browsing was significantly higher for Abies, Fagus, Sorbus and other broadleaf trees species, and lower for Pinus trees in class B (Figure 2). Larix and class A Pinus trees were not significantly different. The explained variance of browsing in model for size class A and in model for size class B was 0.09 and 0.18 respectively.

Effects on Specific Composition
The stress values in the two-dimensional nMDS for all size classes were respectively 0.14, 0.18, 0.16, and 0.16. According to the rule of thumb in Hair et al. (1998) [26], these values indicate that Forests 2020, 11, 612 6 of 12 two dimensions were adequate to represent the variation in measured data. Based on the stress plots, the overall configuration fit well with the data in all the analyses performed (see Figure S3, Supplementary Materials).
For size class A, centroids of the year groups shifted mostly from left to right, along the nMDS1 axis (Figure 3), suggesting that the most important variation in composition across time was the increase in density of European rowan (see Table S2, Supplementary Materials). This variation in species composition between the four years was indeed significant (ANOSIM: R = 0.072, p = 0.001). The nMDS performed on size classes B and C provided similar results. In class B, the shift of the centroids along the nMDS1 axis (Norway spruce = 0.31) suggests that the most important variation in species composition was the decrease of Norway spruce density across the years (ANOSIM: R = 0.01, p = 0.02). Similar changes in specific composition occurred also in class C (ANOSIM: R = 0.015, p = 0.01). In class D, all centroids were located around zero values of both axes and they were close to each other. This suggests that no important shifts occurred in the community during the study (ANOSIM: R = −0.0006, p = 0.49).
explained variance of browsing in model for size class A and in model for size class B was 0.09 and 0.18 respectively.

Effects on Specific Composition
The stress values in the two-dimensional nMDS for all size classes were respectively 0.14, 0.18, 0.16, and 0.16. According to the rule of thumb in Hair et al. (1998) [26], these values indicate that two dimensions were adequate to represent the variation in measured data. Based on the stress plots, the overall configuration fit well with the data in all the analyses performed (see Figure S3, Supplementary Materials).
For size class A, centroids of the year groups shifted mostly from left to right, along the nMDS1 axis (Figure 3), suggesting that the most important variation in composition across time was the increase in density of European rowan (see Table S2, Supplementary Materials). This variation in species composition between the four years was indeed significant (ANOSIM: R = 0.072, p = 0.001). The nMDS performed on size classes B and C provided similar results. In class B, the shift of the centroids along the nMDS1 axis (Norway spruce = 0.31) suggests that the most important variation in species composition was the decrease of Norway spruce density across the years (ANOSIM: R = 0.01, p = 0.02). Similar changes in specific composition occurred also in class C (ANOSIM: R = 0.015, p = 0.01). In class D, all centroids were located around zero values of both axes and they were close to each other. This suggests that no important shifts occurred in the community during the study (ANOSIM: R = -0.0006, p = 0.49).

Damage Analysis
Browsing activity increased from 1994 to 2008 then slowed down in 2014; the overall ungulate abundance shows almost the same tendency [21] (see Figure S2, Supplementary Materials). However, Forests 2020, 11, 612 7 of 12 using the overall abundance of ungulates as the main metric could partially mask the spread of a single species, leading to biased interpretations of the results. The negative trend of chamois population registered between 2003 and 2008 partially compensated for the high increment of red deer. After 2008, the strong depletion of chamois led the overall ungulate population to decrease despite the sustained red deer abundance. These population trends may lead to underestimate the impacts of red deer which has higher nutritional needs compared to the other ungulate species [24,25], and it can use a wider territory under high ungulate density conditions [27]. These inferences, in addition to the partial spatial segregation among chamois and red deer [28], lead us to reasonably assume that most of the browsing marks detected in this study were made by red deer. As a consequence, the reduction of browsing pressure in 2014 may be not linked to lower ungulate abundance, but it could have been caused by other factors, such as extreme meteorological events. The short-lived and thin snowpack registered in winter 2011-2012 (see Figure S4, Supplementary Materials) may have allowed deer to find other food resources, and to reduce the browsing pressure on trees in the regeneration stage [29]. Furthermore, we also hypothesize that the snowfalls that occurred in winter 2013-2014 may have affected deer browsing activity. In that winter, the snow cover at ground level persisted for almost seven months, from November until the end of May. The prolonged snow cover protected the ground vegetation from browsing pressure and, at the same time, it could have reduced the local deer density by leading animals to spend more time outside the reserve or at lower elevations. These hypotheses are both supported by the data collected by weather stations inside the reserve [30]. Due to the lack of empirical evidence concerning the species-specific interaction between a plant and its consumer, we prefer to discuss the results in the following sections referring to ungulates as a guild.
The most browsed size classes were A and B, suggesting the most threatened trees are those between 10 and 150 cm height, as reported in other studies [31]. Taller plants of classes C and D presented lower browsing percentages, indicating that they can be considered mostly out of the ungulate feeding zone and, hence, have more opportunities to reach the canopy.
The clear difference among the species in terms of browsing damage underlines the ability of ungulates to apply different selective behavior on their feeding habits. Apex removal indeed occurred more often on broadleaf species and silver fir compared to Norway spruce, European larch and stone pine. Browsing pressure on selected species was constantly higher than their bearable limits, indicating that many trees would not be able to grow enough to be competitive. Even though Norway spruce was among the least browsed species, it experienced periods of increased browsing pressure. This suggests that this species may act as a resource during food shortage periods [32], or when high population densities force ungulates to feed on lower quality food [33].

Drivers of Browsing Activity
Our investigation about the forces that drive ungulates in feeding on young trees underlined some differences between size class A and size class B plants. Results from the GLM of size class A show that the percentage of adult broadleaf trees in the plot negatively affected the browsing probability. Broadleaf trees in the canopy can decrease browsing pressure on seedlings by reducing the richness of food resources in the ground layer. In contrast, coniferous species in the Pinus and Larix genera promote a diverse understory, then the relative proportion of broadleaved species in mixed deciduous-coniferous stands negatively affects understory species richness [34]. This effect could have led ungulates to spend more time in conifer-rich stands relative to broadleaf-rich stands, especially in summer, when ground vegetation is exposed to herbivory. As a result of this preference, understory vegetation under broadleaf-rich stands is less exposed to browsing pressure. This hypothesis is also confirmed by the results of the size class B GLM, which showed that broadleaf trees did not affect browsing activity, likely because saplings between 30-150 cm of height are browsed in winter, when ground vegetation is covered by the snow [35] and does not attract ungulates. The negative effect of broadleaf trees on browsing pressure of size class A plants may be enhanced by the leaves that drop from deciduous tree sand cover low vegetation, hiding seedlings from ungulates. The effect of aspect suggests that browsing pressure on size class B trees was more intense on northern exposures than on southern exposures. This can be due to the sensitivity of ungulates to snowpack, and to the high ungulate density that may have forced some individuals to forage in areas that are normally avoided. Usually, ungulates prefer to spend wintertime on south-facing slopes, because the climate conditions are milder there, and the thinner snow cover allows them to easily dig and access understory vegetation. The higher availability of food resources may have slightly decreased the browsing pressure on class B trees of these slopes. On the other hand, north-facing slopes are avoided, because the deeper and harder snowpack conditions [36,37] make movement more energy consuming and fewer food resources are available. In these conditions, ungulates like red deer can survive by switching their feeding habits entirely on to class B trees [38], the only food resources available. Moreover, the browsing pressure on north-facing slopes affects class B trees longer than south-facing slopes, because the snowpack is deeper and lasts for more time [39]. We hypothesize that, under high ungulate density conditions, this difference of browsing pressure among different slope orientations may be affected by higher intraspecific competition in preferred areas, which forces some ungulates to frequent northern sites [40].
Both GLMs revealed that browsing pressure was higher in mature forest stands. Browsing probability of size class B plants was positively related to basal area, and negatively associated with tree density, suggesting that ungulates spent more time in mature stands with few large trees than in stands with many smaller trees. These preferences may be due to the difficulties large herbivores face in moving across dense stands, and to greater richness in food resources in the understory layer of mature forests. Indeed, in mature stands, more light reaches the ground, leading to higher seedling and sapling density and higher shrub abundance [41]. This hypothesis is also partially confirmed by the results from GLM for size class A plants, which show that basal area had the same positive effect shown by class B plants, even though tree density did not influence browsing probability.
Lastly, our analysis confirmed that the most powerful driver of browsing pressure was palatability. Species known to be highly palatable such as rowan, beech, silver fir and other broadleaf trees had the highest probability of being browsed. The tender shoots of broadleaf trees and the needles of silver fir are extremely important in the winter diet of ungulates, because they provide energy and nutrients when other food resources are scarce, making them fundamental for survival [42]. On the other hand, species like Norway spruce, stone Pine and European larch have lower digestibility and lower metabolizable energy content, which make them a lower quality food resource [32].
Other factors we did not consider in this study could have played a role in ungulate-forest interactions, such as weather, human disturbance or landscape heterogeneity at larger scales. Fine-scale weather variability (e.g., temperature, moisture, wind chill or snow cover) can strongly affect herbivore movement [43] in ways that are not captured by topographical proxies. Additionally, these effects vary between seasons, i.e., the movement rate of deer is negatively related to cold temperatures during winter, and negatively related to high temperature in summer [44]. Unmeasured human disturbances (e.g., variable hunting pressure, recreation access to forests, or silvicultural activities in logged areas and forest roads) can trigger avoidance responses and force ungulates to change their feeding habits or move to less-disturbed areas [45]. Finally, our sampling design may have failed to capture larger-scale influences of land cover and its heterogeneity, and may be responsible for part of the unexplained variability in our models. An increased interspersion of resources would reduce the distance traveled to gain a particular resource, and, thus, result in smaller levels of movement by deer, especially in winter, when we expect tree browsing to be most intense [46].

Effects on Specific Composition
Since forest regeneration is a long process, the twenty years covered by this study represent a relatively short period of time in which to observe changes. As a consequence of this limitation, we did not expect to observe large differences in terms of changes in species composition. Nevertheless, we were able to detect small variations in the abundance of each species within our plots. Because we did not mark and measure every plant, we cannot determine whether browsing decreases young tree density by directly killing them or by slowing their growth. We suppose that browsing acted in both ways. We hypothesize that seedlings in size class A may die after a single bite, while saplings of size class B are more resistant, but may die later, mainly due to the loss of competitiveness after prolonged exposure to browsing.
The increase in European rowan density in size class A could be caused by the improvement of seed dispersion and the germination processes facilitated by wild ungulates. Indeed, they are known to act as seed dispersal agents by eating fruits, which pass through their guts, and are released on the ground elsewhere (endozoochory) [47]. This dispersal mechanism is operational, particularly in autumn, when many species, such as European rowan, are fruiting. Furthermore, wild ungulates can disperse seeds in other ways, such as by carrying them on their coats or between their hooves (fur-epizoochory and hoof-epizoochory) [48]. On the other hand, in spite of the increasing density of rowan seedlings, we did not find any increase in the frequency of trees belonging to the upper size classes for this species. This suggests that rowan seedlings are not able to grow and establish, due to high browsing pressure. This bottleneck effect causes the failure of recruitment of larger trees, consistent with what has been observed in other studies [13,49].
The dominance of rowan plants in size class A was replaced by the unpalatable Norway spruce in larger size classes. Such dominance of spruce was well maintained over the study period. However, classes B and C showed a small setback in Norway spruce abundance, likely due to the increased browsing pressure registered on classes A and B that constantly exceeded the tolerable limit for this species [16]. The high ungulate population densities during the period of our study probably increased competition for food and forced some animals to become less selective, leading to a higher consumption of unpalatable plant species [50]. No other significant shifts occurred in species composition, suggesting that palatable species had no chance to recover.

Conclusions
This research confirms that overpopulation by ungulates and the resulting heavy browsing pressure on vegetation are powerful forces that shape forest structure and composition. Our results indicate that ungulates browse mostly on trees between 10 and 150 cm of height in mature stands, and that the browsing pressure was consistently higher than the bearable limit for many tree species in our study. We also theorize that browsing on plants in the range of 10-30 cm and browsing on plants between 30-150 cm are seasonally separated events. This study evidences some features of high ungulate density in the study area, such as high amounts of browsing on palatable species, increasing browsing pressure on unpalatable plants, high browsing pressure on less nutritious food resources in unsuitable environments and increasing rowan seeds dispersion. Our results reveal that ungulate overpopulation is indeed currently hampering the regeneration process of European rowan and leading to the dominance of Norway spruce. Moreover, we detected a slight decline in Norway spruce saplings, due to the increased browsing activity of abundant ungulates.
If these impacts continue over the time, they could lead to severe biodiversity losses in the adult tree populations, generating cascading effects on other plant and animal species. Additionally, the reduction of species richness may reduce the resilience of the forest to natural threats, such as fire and windthrows. Further investigations should consider other factors involved in the ungulate-forest interaction. For example, it may be useful to evaluate how the ongoing spread of wolf (Canis lupus L., 1758) will affect the dynamics of browsing, and how it will impact the forest regeneration processes.