The interacting effect of climate change and herbivory can trigger large‐scale transformations of European temperate forests

In many regions of Europe, large wild herbivores alter forest community composition through their foraging preferences, hinder the forest's natural adaptive responses to climate change, and reduce ecosystem resilience. We investigated a widespread European forest type, a mixed forest dominated by Picea abies, which has recently experienced an unprecedented level of disturbance across the continent. Using the forest landscape model iLand, we investigated the combined effect of climate change and herbivory on forest structure, composition, and carbon and identified conditions leading to ecosystem transitions on a 300‐year timescale. Eight climate change scenarios, driven by Representative Concentration Pathways 4.5 and 8.5, combined with three levels of regeneration browsing, were tested. We found that the persistence of the current level of browsing pressure impedes adaptive changes in community composition and sustains the presence of the vulnerable yet less palatable P. abies. These development trajectories were tortuous, characterized by a high disturbance intensity. On the contrary, reduced herbivory initiated a transformation towards the naturally dominant broadleaved species that was associated with an increased forest carbon and a considerably reduced disturbance. The conditions of RCP4.5 combined with high and moderate browsing levels preserved the forest within its reference range of variability, defining the actual boundaries of resilience. The remaining combinations of browsing and climate change led to ecosystem transitions. Under RCP4.5 with browsing effects excluded, the new equilibrium conditions were achieved within 120 years, whereas the stabilization was delayed by 50–100 years under RCP8.5 with higher browsing intensities. We conclude that forests dominated by P. abies are prone to transitions driven by climate change. However, reducing herbivory can set the forest on a stable and predictable trajectory, whereas sustaining the current browsing levels can lead to heightened disturbance activity, extended transition times, and high variability in the target conditions.

Large wild herbivores are an important factor that affects forest demography, community composition, biodiversity (Bernes et al., 2018), and biogeochemical cycles (Leroux et al., 2020).Effects of herbivores are complex, from regeneration browsing and bark stripping, physical disturbance of the soil surface, and seed dispersal, to redistribution of nutrients via urination and defecation (Linnell et al., 2020).Herbivory can trigger cascading effects on tree species composition, ground vegetation, tree regeneration (Bernes et al., 2018;Kuijper, 2011), and different vertebrate and invertebrate communities (Foster et al., 2014).Despite the significance of herbivory increases as game populations become more abundant (Champagne et al., 2021;Morellet et al., 2011), these effects have received considerably less attention compared to stand-replacing disturbances such as windthrows and wildfires.Since the 1960s, there has been a significant increase in herbivore populations in Europe (Apollonio et al., 2017;Spitzer et al., 2020), which mainly concerned Cervus elaphus L., Capreolus capreolus L., Alces alces L., Ovis aries subsp.musimon Pallas, and non-native Cervus nippon Temminck (Carpio et al., 2021;Spitzer et al., 2020).The defining characteristics of overabundance primarily encompass ecological aspects, such as detrimental effects on ecosystems, and socio-economic factors, such as wildlife-human conflicts and damage to crops and forest regeneration (Carpio et al., 2021;Spitzer et al., 2020).The causes of overabundance are complex, including land abandonment and re-naturalization of habitats, restrictive hunting legislation, and the reduction of natural predators (Behr et al., 2017;Spitzer et al., 2020;Valente et al., 2020).Climate change could potentially lead to an increase in ungulate populations by enhancing survival rates and facilitating range expansion (Dawe & Boutin, 2016;Weiskopf et al., 2019); however, these projections remain uncertain due to possible shifts in inter and intraspecific interactions, host-parasite interactions, and other factors (Mason et al., 2014;Morales-Castilla et al., 2021).
Ungulate browsing has been identified as a direct threat to forest adaptation to climate change, particularly by disabling natural adaptive responses of forest vegetation and hindering the achievement of desired species compositions through management (Champagne et al., 2021).Regeneration browsing can shift community composition toward less palatable (generally, with a higher C:N ratio) and browsing-resistant species (Bernard et al., 2017;Fuchs et al., 2021;Katona et al., 2013).Locally, herbivory even causes regeneration failure unless protective measures, such as fencing, are applied (Rhodes et al., 2017).For example, while Quercus spp.and Abies alba Mill.are promoted by forest adaptation strategies in Europe (Hlásny et al., 2014;Lindner et al., 2010), they are also preferred forage for ungulates, which often results in their eradication from species composition (Bernard et al., 2017;Redick & Jacobs, 2020).On the contrary, a widespread Norway spruce (Picea abies [L.] Karst.), which increasingly suffers from climatic stress and biotic disturbance (Hlásny, Zimová, et al., 2021;Zang et al., 2014), is avoided by game.
These effects are context-dependent, affected by tree species identity and diversity, forest age, animals' density and competition, season, and management (Bergqvist et al., 2018;Möst et al., 2015;Tremblay et al., 2007).For example, C. elaphus and C. capreolus in Central Poland favoured European beech (Fagus sylvatica L.) if their populations were high, while the species was avoided at low game densities.In Central Germany, browsing pattern showed a specific scale-dependency, with varying mechanisms driving forage preference and intensity at the regional scale and the patch scale (Ohse et al., 2017).
The increasing pressure from large wild herbivores coincides with the recent disturbance wave in Europe, which culminated after the 2018 drought (Senf & Seidl, 2021).The emerging disturbed areas, often hundreds of hectares large, are exposed to overpopulated ungulates and increasingly hot and dry weather, which collectively affect regeneration vigor and community composition (Vuorinen et al., 2020).This combination of stressors can initiate different recovery pathways and may critically determine future climate regulation, habitat quality, and the provision of ecosystem services (Seidl & Turner, 2022).This particularly concerns P. abies forests, which represent one-quarter of Europe's growing stock.Many of these forests have been intensively cultivated for several centuries outside their native range (Ellenberg et al., 1992), which contributed to their climatic vulnerability and low resilience to natural disturbances (Marini et al., 2012(Marini et al., , 2017)).The recent warm ecological resilience, ecosystem modelling, European temperate forests, forest dynamics, large wildlife herbivores, natural disturbances and dry climate has compromised these forests` defence against insects and pathogens (Netherer et al., 2021), notably the Eurasian spruce bark beetle (Ips typographus L.), and triggered unprecedented eruptions of forest mortality (Hlásny, König, et al., 2021;Patacca et al., 2023).After 2018, ca 4% of spruce growing stock in Europe (8 billion m 3 ) were affected by bark beetles, with the Czech Republic becoming an epicentre.The country has lost 20% of the initial spruce growing stock (500 million m 3 ) after 2018 (Figure 1; Hlásny, Zimová, et al., 2021).
We investigated the disturbance and recovery dynamics of these forests using resilience theory (Holling, 1973) that is increasingly used in forest ecosystem research (Nikinmaa et al., 2020).We interpret resilience as the amount of perturbation (climate change, browsing pressure, etc.) a system can absorb before reaching a threshold beyond which it transitions into an alternative state (Carpenter et al., 2001;Walker et al., 2006).Because critical transitions in forest ecosystems typically occur after years and decades, we employed an approach based on dynamic ecosystem modelling, which is gaining an increasing attention in ecosystem resilience research (Albrich, Rammer, & Seidl, 2020;Egli et al., 2019).
The objectives of this study include determining: (i) the longterm trajectories of forest development under varying degrees of climate change and herbivory, (ii) the corresponding alterations in forest structure, composition, and disturbance dynamics, and (iii) the characteristics of transition pathways triggered by various combinations of stressors and the new equilibrium states towards which the forest converges.We addressed these objectives by the mechanistic landscape-scale model iLand, which has been instrumental in different aspects of forest resilience research (Albrich, Rammer, & Seidl, 2020;Hlásny, Augustynczik, et al., 2021;Rammer et al., 2021).Specifically, we used the model to test the following hypotheses: Climate change combined with an increased disturbance from wind and bark beetle outbreaks will trigger a transition towards climate-adapted forests dominated by broadleaved species, a process that has been envisaged in the scientific literature (Dyderski et al., 2018;Hanewinkel et al., 2013;Thom, Rammer, & Seidl, 2017) and observed locally (Pretzsch & Dieler, 2011).
Alternatively, we hypothesize that the sustained herbivory pres-

| Study area
The study area Kostelec nad Černými lesy (centre coordinates: lat.49°95′, long.14°78′; Figure 2b) is located in the Czech Republic, Central Europe.The size of the area is 40,928 ha with 45% forest cover (18,500 ha).The elevation range is 240-540 m a.s.l., and the annual mean air temperature and precipitation sum ranges are 7.8-9.2°Cand 604-683 mm, respectively, based on 1961-2018 data.The forests have been managed primarily for timber production during the last ca.200 years, which resulted in a widespread occurrence of P. abies, a species native to mountain environments but highly productive across a broad range of conditions (Gömöry et al., 2012).Except for Norway spruce (57.7%), the current tree species composition consists of Scots pine (Pinus sylvestris L.; 10.7%),European beech (7.6%), Sessile oak (Quercus The disturbance rate dramatically increased after 2018 due to the sequence of hot and dry years that triggered the countrywide outbreak of I. typographus that continues until today (Hlásny, Zimová, et al., 2021).Heavy regeneration browsing by overpopulated ungulates is another stressor in place (Vacek, 2017); national estimates indicate two to ten-times higher population levels of C. elaphus, C. capreolus, and Dama dama, and nearly 30-times higher levels of non-native C. nippon, compared to the estimated carrying capacity of the environment (Czech Statistical Office, 2021).The dominant management system is shelterwood (a progressive cutting that leads to the establishment of a new cohort of trees under the canopy of the retained mature trees) with an average rotation length of 120 years.An active search for infested trees, followed by their salvage harvesting, is implemented to prevent the spread of bark beetles.The ownership of the forests is distributed among the Czech University of Life Sciences in Prague (30.1%), the State Forests of the Czech Republic (28.3%), and several hundred private forest owners (41.6%).The landscape includes a beechdominated nature reserve that covers 3.4% of the forest area.

| Simulation model
We used iLand, a process-based forest model that simulates forest demographic processes at the level of trees on landscapes that typically span several thousand hectares (Seidl, Rammer, et al., 2012).Driving environmental variables include temperature, soil water availability, vapor pressure deficit, soil nitrogen availability, and atmospheric CO 2 .A daily soil water balance is computed for a single soil layer accounting for precipitation inputs, interception and evaporation from the canopy, snow deposition and melting, transpiration demand from the canopy, and runoff of excess water.A trees' response to soil water stress is modelled to decrease linearly between the field capacity and a species-specific minimum soil water potential (SWP; Wullschleger & Hanson, 2003).Water availability controls several processes, including photosynthesis, regenerations success, and tree stress mortality (Güneralp & Gertner, 2007).Productivity response to CO 2 is implemented based on the Michaelis-Menten formulation (Friedlingstein et al., 1995).The β factor that expresses the growth response to a relative CO 2 increase was defined based on FACE experiments to 0.6 (Norby et al., 2005).iLand accounts for the amplified fertilization effect under the conditions of water limitation, and downregulation when nitrogen is limiting.However, the nitrogen cycle is not dynamically simulated, thereby not accounting for changes in plant available nitrogen (kg ha −1 year −1 ) due to alterations of vegetation structure, decomposition rates, and atmospheric depositions.Instead, the stable level of nitrogen is determined during the model initialization phase for resource units with a resolution of 100 × 100 m, considering both the input of nitrogen from the soil and the deposition from the atmosphere.The negative effects of nitrogen saturation on vegetation (Tian et al., 2016) are not considered.
Production physiology is modelled using a light-use efficiency approach, driven by environmental conditions and species traits (Landsberg & Waring, 1997).Carbohydrate allocation in trees is calculated annually based on allometric ratios and is sensitive to a tree's competitive status.The mortality probability of a tree is influenced by its carbon balance (stress-related mortality), size, and age.
iLand simulates natural tree regeneration in an annual step in 4 m 2 grid cells.Seed dispersal, climate-dependent establishment, and seedling and sapling growth are driven by species-specific traits (Seidl, Spies, et al., 2012).The dispersal is calculated through speciesspecific dispersal kernels, considering species-specific return intervals of seed years (Lexer & Hönninger, 2001).The establishment is driven by species-specific climatic constraints (minimum winter temperature, minimum and maximum growing degree days, and growing degree days for successful bud burst), soil water potential, and light availability related to emerging canopy openings and modified tree competition.The final establishment is driven by an interaction of seed availability, climatic, and light-related constraints.In case of managed forests, natural regeneration occurs concurrently with tree planting defined based on prescribed species proportions and densities for different sites or arbitrary defined management units (Rammer & Seidl, 2015).
The regeneration can be affected by game browsing implemented through species-specific browsing probabilities (Appendix B in Data S1, Table B1).The probabilities are constant, without considering possible changes in browsing pressure in response to changing vegetation structure, game management, and other factors.
The browsing affects saplings (represented as cohorts rather than individuals) up to the height of 2 m.If the saplings are affected, the cohort does not grow in the given year.If browsing occurs in several consecutive years, the cohort dies.The critical number of years corresponds with the species-specific browsing-tolerance and ranges from a minimum of 2 years, for example, for Sycamore maple (Acer pseudoplatanus L.) and European ash (Fraxinus excelsior L.), to 5 years for European silver fir (Abies alba Mill.).Regarding the most prevalent species in the study landscape, P. abies can endure 5 years of sustained browsing, F. sylvatica 4 years, and P. sylvestris 3 years (see Table B1 for the remaining species).
Management is facilitated through the use of iLand's Agent-Based Management Engine (Rammer & Seidl, 2015), which employs the so-called stand treatment programs (STPs) to execute a series of silvicultural operations throughout the course of stand development.Wind impacts are simulated based on meteorological data such as maximum wind speed, wind direction, and storm duration (Seidl et al., 2014).The disturbance is initiated in locations where vertical differences between the top heights of neighboring grid cells exceed 10 m, i.e., at forest edges (e.g., Blennow & Sallnäs, 2004).A critical wind speed is calculated separately for stem breakage and uprooting based on Gardiner et al. (2000).Wind impacts are simulated iteratively (i.e., wind impacts coming in multiple waves), with forest structure being updated after each iteration.The model was tested, for example, in Sweden following the storm Gudrun, displaying a good correspondence of observed and predicted wind damage patterns (Seidl et al., 2014).

STPs
Bark beetle disturbance considers bark beetle phenology and development, the spatially explicit dispersal of beetles, host tree colonization and defense, and temperature-related overwintering success (Seidl & Rammer, 2017).Non-structural carbohydrate reserves approximate the tree's defense (e.g., Huang et al., 2020).An increase in beetle population size is predominantly triggered by windthrows providing a surplus of breeding material (Wermelinger, 2004)

| Model testing
This is the first study applying iLand in the current study landscape (Figure 1), therefore we thoroughly evaluated the model's performance (Appendix C in Data S1).The test of forest productivity addressed six most abundant species in the study region: P. abies (1514 testing plots), F. sylvatica (337 plots), L. decidua (26 plots), P. sylvestris (313 plots), Q. petraea (123 plots), and A. alba (121 plots).The testing aimed to compare the simulated tree height at a standardized age of 100 years (so-called site index) with index values identified in a field forest inventory (Source: Forest Management Institute of the Czech Republic).The Pearson's correlation of simulated and observed SI ranged from 0.73 (A.alba) to 0.93 (Larix decidua), indicating that the spatial pattern of productivity in the study region was adequately captured.However, the simulated heights were underestimated for all tree species, with an average underestimation of 12%.The mode l's ability to accurately reproduce other indicators of forest productivity (e.g., stand increment, climate sensitivity, tree dimensions) has been demonstrated in previous studies (Albrich et al., 2018;Albrich, Rammer, & Seidl, 2020;Silva Pedro et al., 2015;Thom, Rammer, Dirnböck, et al., 2017).The overall plausibility of simulations was evaluated by comparing a 1000-year long unmanaged and undisturbed simulation with the Natural Potential Vegetation (NPV) of the region that was mapped in the field based on site and phytocenological surveys (Source: Forest Management Institute of the Czech Republic; Appendix C in Data S1, Figure C2).The simulated species composition and the NPV composition mapped in the field (F.sylvatica 47%, Q. petraea 42%, A. alba 7%, other species 4%) were compared using the Euclidean distance similarity, resulting in a 73% match.Finally, the design of simulation plots used for the productivity test was used to evaluate the natural mortality dynamics: the simulated self-thinning coefficients oscillated around the reference values of −1.605 for Reineke and −0.5 for Yoda (Table C2; Reineke, 1933;Yoda, 1963).

| Environmental data
The model requires soil depth, texture, and nitrogen data on a 100 × 100 m grid.These data were derived from the national forest typology map (Source: Forest Management Institute of Czech Republic, n.d.) created based on a field survey.The map contains forest type polygons characterized by a similar site properties and natural potential vegetation, including information about nutrient and hydric conditions, soil depth, structure, and stone content.The available semi-quantitative information was converted into quantitative estimates of effective soil depth and plant-available nitrogen (kg ha −1 year −1 ) according to the internal model logic (Seidl, Rammer, et al., 2012).
Reference climatic conditions were represented by climate data from 1961 to 1990, including daily minimum and maximum temperature, and precipitation.The data were measured at a meteorological station located in the study landscape (lat: 49.9072, long.: 14.7847, elevation: 485 m a.s.l.; source: Czech Hydrometeorological Institute, n.d.).The station measurements were downscaled to the 100 × 100 m environmental grid using the MT-CLIM software (Hungerford et al., 1989).Temperature and precipitation were downscaled using the lapse rates calculated from twelve weather stations in the surrounding of the study region.A vapor pressure deficit was calculated based on temperature and relative humidity data according to Murray (1967).Solar radiation was calculated using MT-CLIM for each grid cell based on its geographical position and topography.Table 1).The climate models were driven by two Representative Concentration Pathways, RCP4.5 and RCP8.5 (Moss et al., 2010).
The four climate models were selected out of 14 models included in FORESEE to represent the variability of projected changes in air TA B L E 1 Projected changes in mean annual temperature, precipitation sum, and vapour pressure deficit for 2071-2100 relative to 1961-1990, according to four climate models driven by two representative concentration pathway scenarios (RCPs).Absolute values for the reference period are indicated in the upper part of the table.temperature and precipitation (Table D1).The daily data from the nearest grid cell (0.1° x 0.1°) of the RCMs were downscaled to each cell of the environmental grid using MT-CLIM.
Since the simulations require 300-year long driving climate data, the reference time series was created by randomly sampling years with replacement from the 1961-1990 period.The atmospheric CO 2 concentration was kept on a level of 332 ppm.For future scenarios, climate data for the period 2021-2100 were defined based on the earlier described climate projections.Data after 2100 were generated by randomly sampling with replacement from the years 2080-2100.The CO 2 concentration up to 2100 was prescribed by the respective RCP scenario (Meinshausen et al., 2011).After 2100, the CO 2 concentration was held constant at the level of 2100, which was 538 ppm for RCP4.5 and 936 ppm for RCP8.5.

| Vegetation initialization
The attributes used to initialize the vegetation in the model were tree species, stand age, and basal area.The data were extracted from

| Simulation experiment
The simulations were run for 300 years, which were found sufficient to evaluate the potential ecosystem transitions.The regional disturbance dynamic was represented by a sequence of windthrows and climatically sensitive outbreaks of spruce bark beetle I. typographus.
The timing of windthrows and wind-related variables (wind speed and direction, windstorm duration; Appendix E in Data S1, Tables E1-E3) were prescribed as a list of events, while outbreaks were simulated as an emergent property of the simulation framework (Section 2.2).
Three 300-year long windstorm sequences were defined to increase simulation robustness.Each sequence consisted of 10 wind events with different properties, defined to reach an average annual disturbance rate of 0.4% (low severity), 0.5% (medium severity) and 0.6% considering the adaptive behavior of forest managers.Details on the management implementation can be found in Appendix F in Data S1.
We defined a set of browsing scenarios and combined them with the climatic scenarios (Table 1, Figure 3).The reference browsing scenario (B reference ) used species-specific browsing probabilities derived from the Austrian forest inventory (originally implemented in iLand) and modified based on different literature sources to better capture the regional browsing patterns (Table B1).Two additional scenarios were created by modifying the reference scenario: a theoretical browsing exclusion scenario (B exclusion ) with all browsing probabilities set to zero, and a moderate browsing scenario with a 50% reduction in reference probabilities (B moderate ), representing successful game control and regeneration protection (Table B1).The total number of simulations was given as: 3 browsing levels × (reference climate + 4 RCMs × 2 RCPs) × 3 wind series, that is, 81 simulation runs.
To assess the effect of climate change and browsing on forest resilience, we adopted the concept of the Historical Range of Variability (HRV), a historical envelope of possible ecosystem conditions under the prevailing disturbance regime (Keane et al., 2009).
This envelope was defined based on the 300-year simulation under reference climate (sampled years from 1961 to 1990), browsing scenario B reference , currently practiced management, and the reference disturbance intensity (Figure E1).Because we focused here on managed forests established ca 200-250 years ago, the thus defined HRV is not fully compatible with HRV typically used in ecological research (e.g., Nagel et al., 2014).Therefore, we use the term reference range of variability (RRV) instead.Forest resilience and transition analysis was based on simulating forest development under all browsing and climate scenarios (Figure 3) and testing them for the departure from the RRV and possible stabilization in alternative states.We conducted this analysis in a two-dimensional ecological space (see e.g., Albrich, Rammer, Turner, et al., 2020;Cushman & McGarigal, 2019;Lamothe et al., 2019) defined by the proportion of conifer species on the landscape and living forest carbon (aboveground and belowground).Considering the proportion of conifers is insightful due to their higher vulnerability to natural disturbances and climate change, which makes it a defining factor of ecosystem resilience in many Central European forests (Honkaniemi et al., 2020;Seidl, Vigl, et al., 2017).The living forest carbon is a sensitive receptor of disturbance impacts and environmental effects on forest growth, productivity, and carbon cycle, making it a frequently used response variable in resilience research (Albrich et al., 2023;Forzieri et al., 2022;Morecroft et al., 2019).The conditions reached at the end of the simulations were referred to as the Target Range of Variability (TRV).
The transition trajectories were characterized by two metrics: the speed and the tortuosity.For this purpose, the two axes of the used ecological space were rescaled into the unit range and the trajectories were smoothed.The speed of trajectory was calculated as the first order derivatives.The point of significant trajectory deceleration, identified based on the second derivatives of the smoothed trajectory, was interpreted as stabilization in new equilibrium conditions.The straightness (tortuosity) was calculated as the ratio of the Euclidian distance between the starting and ending point and the total length of the trajectory; the value of 1 indicates a straight line (Benhamou, 2004).The analyses were conducted using the Rpackage trajr (McLean & Skowron Volponi, 2018).

| Forest structure and composition
Despite the increasingly unfavorable climate, sustaining the reference level of herbivory (B reference ) has maintained a certain proportion of vulnerable yet less palatable P. abies in the species composition (Figure 4b-e; Appendix G in Data S1).Consistent with game foraging preferences (Table B1), drought-tolerant Quercus sp. was removed from the community composition.Meanwhile, F. sylvatica and P. sylvestris increased their proportions and compensated for the loss of P. abies.Forest structure showed remarkable resilience to climate change under B reference : living carbon, mean tree diameter, and tree density (no ha −1 ) did not significantly deviate from the reference conditions under either RCP scenario (Figure 5a,e,h).
On the contrary, reducing herbivory triggered substantial ecosystem reorganization.P. abies declined under both B moderate and B exclusion , while species such as Quercus sp., F. sylvatica P. sylvestris, and Acer pseudoplatanus, which participate in the NPV of the region, gained the dominance (Figure 4c,d,f,g).This transformation was accompanied by substantial changes in forest structure.Living forest carbon decreased under RCP4.5, but it increased under the RCP8.5 by 25%, on average, likely due to the combination of reduced disturbance impact and a strong fertilizing effect of elevated CO 2 concentration (Figure 5b,c).Forest structure was shifted towards smaller-diameter trees with a higher stand density, which mainly resulted from the emergence of a cohort of advanced regeneration (Figure 5f,g,I,j).

| Disturbance dynamics
Although the reference herbivory level (B reference ) helped to preserve the main characteristics of forest structure and composition under climate change (Figures 4 and 5), this development was characterized by a substantial increase in disturbance (Figure 6b).The cumulative disturbance impact was 9% and 40% higher under RCP4.

| Forest resilience and transition
Considering the ecological space defined by the proportion of conifers and living forest carbon, the reference browsing scenario (B reference ) combined with the reference climate depicted a welldefined Reference Range of Variability (RRV), ranging between 68%-80% of conifers and 90-155 tC ha −1 (Figure 7a).Four out of six tested combinations of browsing and climate change departed from the RRV and transitioned to new equilibrium conditions: all browsing scenarios combined with the RCP8.5, and B exclusion combined with RCP4.5 (Figure 7g).Under B reference and B moderate , the Target Range of Variability was considerably larger than under the RRV, suggesting instability of future conditions.This was mainly due to the higher proportion of conifers, which fluctuated over time and across climate models from 15% to 50%.
The browsing exclusion scenario showed the distinct decelera-

| Community composition and forest structure
The investigated interactions between climate change and forest regeneration browsing produced both compounding and antagonistic effects on forest community composition, tree size structure, and carbon.A reference development that maintained the current level of herbivory and considered climate from 1961 to 1990 preserved the main forest attributes close to the initial conditions from 2014, aligning with more than 200-year existence of these forests (Ellenberg et al., 1992;Johann et al., 2004).Exposing the forest to climate change in the simulations, while keeping the reference browsing intensity, led to a reduction in the proportion P. abies.
This aligns with a high climatic sensitivity of the complex P. abies -I.typographus (Marini et al., 2012(Marini et al., , 2017)), recent observations, and model projections indicating a decline of European spruce forests under climate change (Dyderski et al., 2018;Hanewinkel et al., 2013;Hlásny, Zimová, et al., 2021).The reduction of P. abies was compensated mainly by an increase in F. sylvatica and P. sylvestris, species that exhibit higher drought tolerance than P. abies (Pretzsch et al., 2020), are less vulnerable to wind (Wallentin & Nilsson, 2014), ) (a) Klopčič et al., 2010).Notably, oaks, which constitute an important part of the landscapes' NPV, and are favoured by climate change (Mette et al., 2013), were excluded from the emergent species pool due to their high palatability (Petersson et al., 2019).This indicates that selective browsing pressure could hinder the expected adaptive shifts in community composition in response to climate change, as it has already been observed in some regions in North America and Europe (Beguin et al., 2016;Champagne et al., 2021;Côté et al., 2004).However, none of the tested combinations of browsing pressure with climate change resulted in the initially hypothesized critical loss of resilience characterized by a persisting reduction in forest carbon or a regeneration failure, as it has been identified, for example, in some fire-prone systems (Enright et al., 2015;Hansen et al., 2018).
Relaxing browsing pressure (B moderate and B exclusion ), combined with both RCP4.5 and RCP8.5, unlocked adaptive changes in forest structure and composition, resulting in the substantial increase in naturally dominating Quercus spp.A decline in P. abies was more pronounced than under B reference , nearly disappearing from the species composition within 100-120 years (Figure 4c,d,f,g).Such changes in community composition were also identified through a simulation experiment from Germany; however, the overall decline in forest basal area caused by climate change was not compensated by the changing community composition due to the browsing reduction (Cailleret et al., 2014).At the same time, forest structure underwent adaptive changes, being shifted towards a larger number of individuals with smaller diameters.These directions of change agree with previous studies suggesting the expansion of oak-dominated forest types across Central Europe (Dyderski et al., 2018;Hanewinkel et al., 2013), a replacement of P. abies by F. sylvatica (Pretzsch & Dieler, 2011), and a global shift of forest structure towards smaller sizes (McDowell et al., 2020;Stovall et al., 2019).
Under all browsing scenarios, forest response to RCP8.5 was characterized by an increase in living carbon.Although we did not F I G U R E 5 Relative changes in living carbon, mean tree diameter, and the number of trees with a height greater than 4 m, compared to reference simulations driven by a combination of reference browsing intensity and reference climatic conditions.The mean and minimummaximum ranges of simulations driven by four climate models and three replicates of wind disturbance series are shown.The reference values indicate a 300-year average of the simulation.ining the sensitivity of forest carbon recovery to various CO 2 levels using iLand supports such interpretation (Dobor et al., 2018).The identified increase in carbon aligns with the recent global trend of enhanced forest growth attributed to CO 2 fertilization and an extended growing season (Pretzsch et al., 2014;Yuan et al., 2019).
However, the continued strength of this effect is uncertain, poten- fertilization, global plant biomass could see an increase of 12% by 2100.This is lower than the average increase of ca 25% reported in our current study (Figure 5c).An explanation could be a simplified representation of nitrogen cycle in iLand (Seidl, Rammer, et al., 2012), whereas nitrogen availability is crucial for sustaining the fertilization effect (Schimel et al., 2015;Terrer et al., 2019) and significantly affects the outcomes of terrestrial carbon cycle models (Piao et al., 2013).Moreover, the magnitude of growth response to a relative increase in CO 2 was set based on Norby et al. (2005), while some recent studies that more comprehensively accounted for fertilization constraints and plant acclimation (e.g., Keenan et al., 2016), were more conservative in this regard.Other model limitations, such as the lack of consideration for the adverse effects of nitrogen saturation (Tian et al., 2016) and the legacy mortality after droughts due to cavitation (Anderegg, 2015) may have also influenced the results.
These limitations, along with the significant uncertainty of future CO 2 fertilization (Smith et al., 2016), necessitate caution and underscore the need for a more comprehensive integration of ecosystem experiments into vegetation models (Bonan & Doney, 2018;Medlyn et al., 2015).
Finally, while we have considered the intensification of I. typographus outbreaks due to climate change, this effect diminished as P.
abies (the host of I. typographus) declined on the landscape.However, we did not introduce any biotic agent that could impact the emerging pool species, such as Quercus spp., F. sylvatica, and P. sylvestris, neglecting the potential future shifts in biotic disturbance regimes (Forzieri et al., 2021;Seidl, Vigl, et al., 2017) and an accelerating rate of invasions (Seidl et al., 2018).This assumption could also have contributed to the simulated increase in forest carbon.

| Disturbance dynamics
The simulated disturbance impacts distinctly varied depending on the degree of climate change and browsing but also exhibited a specific temporal progression, responding to the transient changes in forest structure, community composition, and trait assemblages.
Under the reference browsing scenario, the disturbance impact was significantly amplified by climate change (Figure 6b), mainly accounting for the warming-induced increase in the number of bark beetle generations, decrease in winter mortality, and the compromised tree resistance (Kausrud et al., 2012).This development aligns with the recent surge in bark beetle outbreaks across Europe triggered by a series of warm and dry years (Hlásny, Zimová, et al., 2021;Senf & Seidl, 2021).A relatively high level of disturbance was sus-

| Societal and ecosystem management implications
Managing ecosystems in transition requires sustained efforts from managers (Millar et al., 2007), including determining when to bolster resilience, when to accept the transitions, and when to actively steer their direction and pace (Schuurman et al., 2022) FCCC, 2015).On the downside, shifting forest demography towards smaller diameter trees and changing species composition from conifers to broadleaves can reduce timber production and the amount of habitat related to the presence of large coniferous trees (Brèteau-Amores et al., 2023;Felton et al., 2016).
Still, although the reduction of browsing can set the forest on an adaptation trajectory, the continued overabundance of game remains a realistic possibility due to logistical, economic, legislative, and societal barriers to game regulation (Carpio et al., 2021).This development may result in high severity of disturbance and unpredictable development trajectories, challenging the limited societal capacity to deal with supply and demand shocks (Millar et al., 2007;Morland et al., 2018).These facts highlight the trade-off between managing forests for resilience, biodiversity, and other ecosystem services, which should be increasingly considered in management and policies (Newton, 2016).

| Methodological considerations
While forest resilience can be significantly shaped by management (Messier et al., 2013), our study has only incorporated a single management scenario.This scenario, while reflecting local practices, also represents a management approach commonly applied across extensive forest areas in temperate Europe (Mason et al., 2022).While such management implementation enhances the relevance of our findings, it neglects the possible interaction of browsing manipulation with alternative silvicultural strategies, which could decisively affect ecosystem resilience (Messier et al., 2019).Furthermore, although the timing of management operations within the utilized STPs has been dynamically adjusted to changing forest conditions, the structure of STPs (Appendix F in Data S1) remained unchanged throughout the simulations.This strategy allowed us to isolate the combined effect of browsing and climate change, but it did not account for societal adaptability to environmental changes through management modifications (Nikinmaa et al., 2023).Therefore, expanding our current research to consider alternative and adaptive management practices could be an important progression towards understanding the socio-ecological dimension of resilience of European temperate forests (see, for example, Nikinmaa et al., 2020).
Defining reference conditions is an essential component of resilience assessment (Lamothe et al., 2019), yet it is challenging in managed forests, which have often replaced previous natural forests in their recent history.This hamper reconstructing the historical disturbance regime (a central component of Historical Range of Variability; Keane et al., 2009) that can be manifested on time scales longer than the investigated forest type exists.For example, the return interval of a large stand-replacing disturbance in Central European mountain forests is more than 300 years (Čada et al., 2020), while here investigated forests were established 200-250 years ago.Therefore, the RRV employed here should be viewed as a theoretical, yet practical construct rather than a rigorous model derived from, for example, dendrochronological reconstructions (e.g., Čada et al., 2013).Conceptually, the RRV defined here is similar to the Simulated Range of Variability (SSR; Cushman & McGarigal, 2019), which highlights that the range is subject to the limitations of the model.
Finally, although the used browsing coefficients corresponded with the reviewed literature, they neglected the frequently reported context-specificity of browsing intensity that depends on the diversity and species composition of the regeneration cohort (Ohse et al., 2017), density-dependent foraging preferences (Bergqvist et al., 2018;Tremblay et al., 2007), and forest and wildlife management (Hothorn & Müller, 2010;Möst et al., 2015).Moreover, the relaxation of temperature and snow cover limitations due to climate change can have a positive effect on ungulates (Rivrud et al., 2019), while high warming levels combined with substantial vegetation transformations can induce behavioural shifts, potentially changing herbivory patterns.Therefore, coupling vegetation models with large herbivore population models (e.g., Owen-Smith, 2009), including effects of predators such as wolves and lynxes that are reestablishing in many parts of Europe, seems an important research avenue in ecosystem modelling, helping resolve some questions of ecosystem resilience and transition.
sure will compromise such a natural adaptation by maintaining the dominance of less palatable yet vulnerable P. abies.Exposure to F I G U R E 1 Areas affected by the large-scale bark beetle outbreak in the Czech Republic after 2018 (a, b); a herd of fallow deer Dama dama (c, Photo: J. Červený); and a replanted site protected from ungulates by a fence (d, Photo: R. Modlinger).climate change could then trigger high disturbance activity, locking the forest in a cycle of disturbance impacts and disabled adaptive responses.
The model can simulate different management operations and disturbances such as wind, bark beetles, and fire in a spatially explicit processed-based manner.iLand uses several temporal and spatial scales: (i) resource availability is simulated for grid cells with a resolution of 100 × 100 m; (ii) resource competition is simulated at the level of individual trees; (iii) environmental constraints are considered on F I G U R E 2 Study area Kostelec nad Černými lesy: (a) location in Europe (red circle) with forest cover in the background (Source: CORINE Land Cover, EEA, 2018).Country codes according to the European norm ISO-3166 Alpha-2 are displayed; (b) study area's elevation and forest cover.Map lines delineate study areas and do not necessarily depict accepted national boundaries.a daily scale; (iv) primary production and soil processes operate on a monthly scale; (v) and mortality, disturbances, and management are considered annually.
include planting, thinning, harvesting, and post-disturbance salvaging.Different timing and intensities of these operations can be defined to address different sites conditions, species composition, and management objectives.The landscape-level scheduling is applied to reach a more balanced interannual harvest distribution and account for spatial contingencies.The adaptive elements of STPs include modifying the timing of thinning based on the actual stand density, replanting disturbed sites depending on the state of natural regeneration, and resetting the sequence of operations if a stand experiences disturbance.
. Minor fluctuations in population size occur irrespective of windthrows depending on climate-sensitive background outbreak probability.Simulated infestations were tested against independent data by Sommerfeld et al. (2021), who found a good correspondence between the simulated levels of infested growing stock and spatial infestation patterns.
Future climate was represented by eight climate projections produced by four regional climate models (RCMs) from the FORESEE 4.0 database (FORESEE, n.d.;Kern et al., 2024).The dataset provides bias corrected climate projections for Central Europe based on the projections from the CORDEX initiative(Jacob et al., 2014; the regional Forest Management Plans (FMP), which were updated in the field between 2010 and 2014.The data represent the average characteristics of forest stands, which are polygons with a mean size 1.04 ha.The data were available for the state-owned forests only, which represent 58.4% of the forest area (13,700 forest stands).For private forests, we approximated stand data based on the categories 'coniferous', 'broadleaved', and 'mixed' from the seamless Europewide landcover classification CORINE Land Cover(EEA, 2018).The landcover polygons were divided into 100 × 100 m cells aligned with the environmental grid.Each cell was populated with tree species, stand age, and basal area data from adjacent forest stands with available stand data (i.e., from the state-owned forests), considering the CORINE Land Cover class, soil, and elevation constraints.Rather than using the FMP data directly for initializing the vegetation, we used a 1500-year long spin-up simulation to ensure that vegetation and environmental conditions were in equilibrium at the beginning of the simulations.The spin-up started from the bare ground with unlimited seed availability and site conditions described in the previous section.The driving climate data were sampled from the observations from the period 1961-1990.The spin-up procedure assimilated information about each stand's age, basal area, and species composition from the FMP database.The spin-up was stopped individually for each stand after converging to the target values, provided the spin-up duration was 300 years, at least.Such a procedure allowed for reaching the initial forest state consistent with the model logic and reliably representing the current structure and composition of the forest(Thom et al., 2018).

(
high severity) of the growing stock.These values characterized the observed disturbance level in 1960-2000.The baseline disturbance definitions corresponded with the reference climate, management, and browsing patterns (see below).Simulations driven by different browsing and climatic scenarios (Figure 3) were exposed to the same wind sequences, but the resulting disturbance impact depended on actual forest structure and composition.The currently practiced management (a shelterwood cut with a rotation length 80-140 years) was introduced based on five stand treatment programs (STPs;Rammer & Seidl, 2015).The STPs were defined for the prevalent site types in the study area (Appendix F in Data S1) in a way that allowed for their implementation irrespective of potential shifts in species composition during the simulation, not F I G U R E 3 Simulation design with three levels of forest regeneration browsing, eight climatic projections (four Regional Climate Model projections driven by two Representative Concentration Pathways, RCPs), a reference climate, and three replicated wind series.
5 and RCP8.5, respectively, than under reference climatic conditions.The largest singular impacts (FiguresE2-E4) affected the growing stock 1.4 and 2.5 times more under RCP4.5 and RCP8.5, respectively, than under the reference climate.On the contrary, changes in forest structure and composition initiated by reduced herbivory (i.e., B moderate and B exclusion ) mitigated disturbance impacts and reduced the climatic sensitivity of the disturbance regime.Specifically, the cumulative disturbance impact was 52%-53% lower under B moderate and 57% and 66% lower under B exclusion compared to the reference simulation (Figure6c,d).At the same time, the difference between RCP4.5 and RCP8.5 was considerably reduced.Finally, the time frame with the most significant disturbance impacts, as indicated by the steepness of the cumulative amounts of disturbed growing stock (Figure6b-d), was notably shortened as the browsing pressure was reduced.
of trajectory speed after 120 and 170 years under RCP4.5 and RCP8.5, respectively, indicating the time to stabilization in new equilibrium conditions (Figure 8d-g).At the same time, the trajectory under RCP8.5 was the straightest (as indicated by the Straightness Index 0.52-0.80, Figure 7g), corresponding with the lack of erratic fluctuations due to the reduced disturbance impact (Figure 6d).The time to stabilization was longer under RCP8.5 with both B reference (215 years) and B moderate (180 years), and the trajectories were more tortuous (Figure 8e,f).No distinct point of deceleration was identified in the remaining simulations (Figure 8b,c).
are only moderately favored by game (Heuzé et al., 2005; F I G U R E 4 Simulated tree species composition, in terms of the living carbon stock, under different levels of browsing and climate change.The reference simulation was driven by a combination of climate corresponding to the period 1961-1990, actual browsing intensity, and the medium severity wind scenario.The effect of climate change was demonstrated based on the simulation driven by a single RCM projection (EC-EARTH_RACMO22E_r1; Appendix D in Data S1) with two RCP scenarios.Simulations driven by the remaining climate models are shown in Appendix G in Data S1, Figures G1-G3.
conduct an effect attribution experiment, this increase was obviously facilitated by the fertilizing effect of CO 2 that compensated for the adverse effects of climate change.A previous analysis exam- tially being compromised due to the escalating counter-effects of increasing VPD and forest mortality(Peñuelas et al., 2017;Yuan et al., 2019).Still,Terrer et al. (2019) proposed that due to CO 2 F I G U R E 6 Cumulative amount of forest growing stock affected by the disturbance under different browsing and climatic scenarios during the 300-year simulation.Grey rectangles indicate the period with largest disturbance impacts.The horizontal dashed line indicates the maximum of the reference simulation (reference browsing × reference climate; panel a).The envelope in panel (a) is based on three wind scenarios.The envelopes in panels (b), (c), (d) are based on four climate models and three wind scenarios.Trajectories of forest development over a 300-year period under different intensities of forest regeneration browsing (columns) and climate change (rows) (b-g).The red dot indicates the starting point of all simulations (forest conditions from 2010 to 2014).The dashed black circle represents the reference range of variability (RRV), defined as the envelope of forest conditions driven by the reference level of browsing and reference climate (a).The red ellipses represent the envelope of target states after 300 years of forest development.The interquartile range of the Straightness Index is shown, with a value of 1 indicating a straight line.
tained during the simulation due to the persisting presence of P. abies, a result of the suppressed competition from non-host species caused by browsing.From a temporal viewpoint, the timeframe between simulation years 80 and 180 experienced the most severe disturbance impacts (Figure6b).This severity came from the strong amplification of the disturbance regime by climate change and the lack of adaptive reorganization within the forest.Later, forest communities shifted towards structures and trait assemblages better suited to the new disturbance regimes (such as smaller statures,F I G U R E 8The speed of trajectories evaluated in a two-dimensional ecological space defined by the proportion of conifers and living forest carbon.A distinct decrease in speed (a vertical red line) indicates trajectory stabilization in new equilibrium conditions.Mean and minimum-maximum envelopes based on three replicates of driving wind series (a), and three replicates of wind series and four climate models (b-g) are shown., etc.), conferring increased resistance to disturbances(McDowell et al., 2020).However, given the ongoing wave of disturbances in Europe(Patacca et al., 2023), this adaptive cycle may be much faster than our simulations have indicated, with disturbances functioning as catalysts of forest adaptation to climate change and new disturbance regimes(Thom, Rammer, Dirnböck, et al., 2017).Finally, forest transformations facilitated by a reduction in browsing have effectively counteracted the amplification of disturbance regime caused by climate change and accelerated the described adaptive cycle.Reducing and excluding the effect of browsing also shortened the period of disturbance culmination (Figure 6c,d) and decreased the variability of disturbance impacts, steering the forests towards a more stable trajectory compared to maintaining the reference browsing level.4.3 | Forest resilience and transitionUnder RCP4.5, the transitions were triggered by the exclusion of browsing only, while the transitions occurred under RCP8.5, regardless of the intensity of browsing.At the same time, high and moderate browsing pressures (B reference and B moderate ) have prevented the transition under RCP4.5, maintaining the forest within its RRV.This combination of stressors thus establishes the actual boundaries of forest resilience, provided there are no substantial alterations in management and disturbance regimes.Yet, this comes at a cost of maintaining the forest in a maladapted state (sensuMalhi et al., 2020), increasing the likelihood of large disturbance impacts in the future.These findings imply that the studied forests are resilient to past disturbance regimes and moderate climate change, which is likely the case for most European forests(Senf & Seidl, 2021).However, our findings also highlight the risks of resilience erosion and substantial shifts in equilibrium forest composition and size structure under aggravating climate change.Furthermore, we revealed that the time required to stabilize in new equilibrium conditions was positively associated with the level of combined stress from browsing and climate change.A relatively minor level of stress, corresponding with B exclusion and RCP4.5, led to the stabilization within 120 years.However, under RCP8.5, the transition times increased from 170 years under B exclusion to 180 years under B moderate to 220 years under B reference .Collectively, the identified boundaries of ecosystem resilience, along with the revealed association of stress levels with transition time, provide a novel perspective on the resilience of European temperate forests.
. Our research has indicated that the studied forests are prone to transition and there are compelling reasons to actively guide them.Apart from stabilizing forest trajectories and reducing disturbance impacts, the transitions initiated by browsing reduction, including increase in forest carbon stocks, align with European climate change adaptation strategies (European Commission, 2021; Lindner et al., 2010) and the goals of the Paris Agreement, which promotes forest management as a means of mitigating climate warming by reducing CO 2 concentrations (UN