Regeneration dynamics in mixed mountain forests at their natural geographical distribution range in the Western Rhodopes

barely known, and their dependence on microsite and management effects needs to be better understood. The objective of this study was to investigate ecological factors under management regimes of different intensity ( “ single-tree ” -selection and “ group-tree ” -selection) that influence the regeneration processes in mixed mountain forests in the Bulgarian Rhodopes. Data on regeneration and microsite conditions were collected on 105 systematically distributed plots (25 m 2 /100 m 2 ) in four 100 – 120 years old stands located in the regional forest district of Smoljan, Bulgaria (1580 – 1650 m a.s.l.). We relied on generalizeds linear mixed models to analyse for each species the (1) size-dependent regeneration density and (2) height increment in dependence on management practices, competing vegetation, as well as soil and light conditions. Our study revealed an overall high potential for recruitment in the Western Rhodopes. Regeneration density was highest in fir (median 12800 N ha-1), followed by spruce (median 1600 N ha-1) and beech (median 1200 N ha-1). Fir benefited most from “ single-tree ” selection cuttings, while “ group-tree ” selection cutting tended to promote beech and fir but also spruce. Competing ground vegetation was detrimental for seedling density of all species. Annual height increment increased with plant size, was lowest in spruce, and similar in fir and beech. Sapling increment was driven by light, whereas seedlings did not react to increased radiation. Browsing was species-specific and was highest in beech (15 – 30 %), followed by fir (5 – 10 %) and spruce ( < 1 %). It was not a crucial factor in impeding tree recruitment. We conclude that frequent harvest activities of low intensity which consider advanced regeneration are a promising approach to successfully convert the formerly spruce-fir-dominated forests to climate-adapted fir-beech-(spruce)-mixed stands.

Mixed mountain forests consisting of Norway spruce (Picea abies (L.) Karst.),European beech (Fagus sylvatica L.), and silver fir (Abies alba Mill.) are among the most productive and stable forest ecosystems in Europe.Their southeasternmost geographical distribution range is located in the Western Rhodopes, where they have high economic, recreational, and ecological value.In the past, shelterwood cuttings dominated forest management practices in these forests and were mainly aimed at maintaining and reproducing conifers.During the past two decades, single-tree and group-tree selection systems have been promoted as alternative management approaches to support the conversion of spruce-dominated stands to close-to-nature mixed forests of fir, beech, and spruce.However, the natural regeneration dynamics in these stands are barely known, and their dependence on microsite and management effects needs to be better understood.
The objective of this study was to investigate ecological factors under management regimes of different intensity ("single-tree"-selection and "group-tree"-selection) that influence the regeneration processes in mixed mountain forests in the Bulgarian Rhodopes.Data on regeneration and microsite conditions were collected on 105 systematically distributed plots (25 m 2 /100 m 2 ) in four 100-120 years old stands located in the regional forest district of Smoljan, Bulgaria (1580-1650 m a.s.l.).We relied on generalizeds linear mixed models to analyse for each species the (1) size-dependent regeneration density and (2) height increment in dependence on management practices, competing vegetation, as well as soil and light conditions.
Our study revealed an overall high potential for recruitment in the Western Rhodopes.Regeneration density was highest in fir (median 12800 N ha-1), followed by spruce (median 1600 N ha-1) and beech (median 1200 N ha-1).Fir benefited most from "single-tree" selection cuttings, while "group-tree" selection cutting tended to promote beech and fir but also spruce.Competing ground vegetation was detrimental for seedling density of all species.Annual height increment increased with plant size, was lowest in spruce, and similar in fir and beech.Sapling increment was driven by light, whereas seedlings did not react to increased radiation.Browsing was species-specific and was highest in beech (15-30 %), followed by fir (5-10 %) and spruce (<1 %).It was not a crucial factor in impeding tree recruitment.We conclude that frequent harvest activities of low intensity which consider advanced regeneration are a promising approach to successfully convert the formerly spruce-firdominated forests to climate-adapted fir-beech-(spruce)-mixed stands.

Introduction
In the past decade, research on ecology and management of mixedspecies forests has proliferated (Albrich et al., 2020;Gillerot et al., 2021;Pretzsch et al., 2017), as it is essential to understand how to adapt forests to climate change impacts (Pluess et al., 2016).In addition, mixed forests are expected to provide a greater range of ecosystem goods and services compared to monospecific forest stands (Bauhus et al., 2017;Gamfeldt et al., 2013).In the temperate zone of Europe, mixed mountain forests consisting of Norway spruce (Picea abies [L.] H. Karst; "spruce"), silver fir (Abies alba Mill.; "fir"), and European beech (Fagus sylvatica L.; "beech") are widely distributed and cover an area of more than 10 million hectares (Hilmers et al., 2019).Europe's forests are strongly shaped by centuries of human activities (Conedera et al., 2017;Kaplan et al., 2009), and most of the mixed mountain forests accessible to humans were often intensively exploited and/or converted to monospecific age-class systems dominated by spruce (Pretzsch et al., 2020).The observed climatic changes since the second half of the 20th century (IPCC, 2021) and the associated increasing uncertainties in the future forest development and health necessitated a re-assessment of such silvicultural practices (e.g.Bowditch et al., 2020;Knoke et al., 2021;Rigling et al., 2008;Seidl et al., 2011).This is especially true for forests dominated by spruce, a species considered to be particularly vulnerable to increasing drought stress (Nikolova et al., 2011;Schmied et al., 2023) and large-scale disturbances that are already causing excessive mortality (Bigler et al., 2006;Gimmi et al., 2010;Senf et al., 2020).Therefore, the conversion of spruce-dominated stands to close-tonature mixed conifer and broad-leaved forests is crucial for forest practitioners (Brang et al., 2014;Hlásny et al., 2017;Pretzsch and Zenner, 2017) and is therefore focused by forest policy makers (BMEL, 2020) and research (Hilmers et al., 2020).
Silvicultural concepts for mixed stands usually aim to establish structural diversity and broad species composition (Pretzsch et al., 2021;Schütz, 1999) by applying single-tree and group-tree selection approaches (Brang et al., 2014).Some knowledge on the transformation process of conifer-dominated forests is now available, mostly from Central European forests (e.g.Hlásny et al., 2017;von Lüpke et al., 2004).However, the existing guidelines for practical use are still largely based on theoretical schemes (Hilmers et al., 2020) and need to be more widely supported by data from other regions and practical examples for conversion (Pretzsch, 2020;Reventlow et al., 2021).
The effectiveness of conversion is coupled with the regeneration process (Axer et al., 2022;Löf et al., 2023), as it determines the adaptive capacity of forests to environmental change and affects ecosystem services supply (Aquilué et al., 2021;Duveneck and Scheller, 2015;Thom et al., 2022).Understanding how light and other microsite conditions influence the species-specific regeneration density and growth dynamics is essential for the recommendation of appropriate management strategies that ensure long-term structured and mixed forests (de Frenne et al., 2021).The ecological requirements of trees differ depending on species and developmental stage (Diaci et al., 2020;Unkule et al., 2022).Height growth is a key determinant of tree recruitment (i.e.development of young trees across a certain size threshold) and is driven by light conditions and species-specific shade-tolerance (e.g.Käber et al., 2021;Klopčič et al., 2015).Light environment consists of direct and indirect solar radiation, both of which contribute to photosynthesis via their photosynthetically active photon flux density (PPFD) (Grasovsky, 1929;Muraoka et al., 2001).However, they may have different effects on factors determining microsite quality, such as temperature or waterbalance, which is particularly relevant during seedling establishment (e.g., Brang, 1998;Diaci et al., 2020;Imbeck and Ott, 1987).Regeneration dynamics is considered to reflect species-specific abilities of plants to transform favourable light conditions into growth and their survivorship under the canopy.This is described as shade-tolerance (Kobe et al., 1995) and can change, species-specific, during tree ontogeny (Kneeshaw et al., 2006;Petrovska et al., 2022).Therefore, understanding the role of factors controlling natural regeneration processes in mixed mountain forests is a field of intense research (e.g.Ammer and Weber, 1999;Bauhus et al., 2017;Löf et al., 2018).While most of the mixed spruce-fir-beech forests are naturally located in the Alpine arc, they are also common in mountainous regions of Eastern Europe (Hilmers et al., 2019;Horvat et al., 1974;Leuschner et al., 2017).The south-easternmost distribution limit of this species composition is in the Western Rhodopes, a mountain range located in southern Bulgaria.At an altitude between 1000 and 1600 m a.s.l., climatic conditions are optimal for beech and fir, while spruce reaches its warm temperature limit (Kölling, 2007;Zlatanov et al., 2017).In addition, weather extremes such as dry summer droughts and heavy rainfall events are already occurring more intensively and frequently in the Balkans (Aleksandrov et al., 2009;Bocheva et al., 2009).These extremes likely further weaken the stability of some spruce-dominated forests (Zlatanov et al., 2017).Until the 1990s, forest management in the Western Rhodopes was characterized by shelterwood cuttings, with spruce being intentionally prioritized for more valuable timber and for recreational purposes (Zlatanov et al., 2017).In the past two decades, however, principles of "close-to-nature" silviculture (Schütz, 1999) have been applied more frequently.Thus, forest practices have shifted from traditional production-oriented forest management to alternative silvicultural approaches to promote stand heterogeneity.As the mountain forests in this region are an essential economic factor for the local population and are crucial for the recreation, biodiversity, erosion control, and water quality of the region (Panayotov et al., 2019;Zlatanov et al., 2017), the successful conversion of formerly sprucedominated stands to mixed stands of conifers and broad-leaved trees is aspired by foresters.
To our knowledge, however, there are no studies focusing on the natural regeneration processes during early stages of conversion of spruce-dominated stands, in particular at their south-eastern distribution margins.To address this gap in knowledge, we conducted our study in four formerly spruce-dominated stands: two of which were recently treated by moderate thinning ("single-tree"-selection) and the other two by heavy thinning ("group-tree"-selection) for developing structural heterogeneity more rapidly.The four study stands were located on north-faced slopes in the Western Rhodopes, were of similar age, and developed under comparable soil and climatic conditions.Our main questions were: 1) What are the main factors driving the stem density of fir, spruce, and beech regeneration in formerly spruce-dominated stands undergoing conversion to mixed fir-spruce-beech forests?
2) To what extent is the height growth of young trees dependent on light conditions and how does this process vary between tree species and sapling size?

Description of the study area
The study area is located in southern Bulgaria in the Western Rhodopes.This mountain range covers an area of about 12200 km 2 , of which 8700 km 2 are on elevations between 1000 m a.s.l. and 2191 m a.s.l (Heiss and Josifov, 1990;Shishkov and Kolev, 2014).
In this geographical region, spruce grows naturally in a mixture with beech and fir, mostly on aspects with northern exposure (Panayotov et al., 2016).Such mixed forests are located at the south-easternmost edge of their natural range (Fig. 1 and Fig. S1).The climate is montane (Dimitrov, 1976;Grunewald et al., 2009;van Huis et al., 2013) with mild winters, reflecting the Mediterranean influence and humid summer months (Fig. 1).The mean annual temperature is 8.8 • C and the annual precipitation sum is 920 mm, as shown in the long-term records of the climate Station Chepelare (Fig. 1) for the period 1990-2019 (Fig. 1).In the context of the local climate for the period 1940-2019, the last 30 years were characterized by increasing air temperatures with D. Ambs et al.
The soils of the region are described as brown forest soils, or Cambisols, which developed over rhyolite, and provide sufficient water and nutrient supply (IUSS Working Group WRB, 2015;Panagos et al., 2011;Shishkov and Kolev, 2014).According to the forest management plans, the four study sites have similar soil and climatic characteristics.

Selection and characteristics of the study sites
The field work was conducted in two forestry districts Pamporovo and Smolyan (Fig. 1) in forests with high recreational and economic value.Within each forest district, two study stands were chosen: PAM1 and PAM2 in Pamporovo, and SMO1 and SMO2 in Smoljan.All study stands were selected according to three criteria: (1) dominating conifers in the overstory (basal area share of spruce and fir ≥70 % of total basal area), (2) age of the trees in the overstory >100 years, and (3) period since the last harvest activity ≥10 years.All stands were pre-selected based on forest management plans from the year 2019 and then verified by field observations.
The study stands were located at elevations of 1580-1650 m a.s.l. on north-west to north-east facing slopes and have an inclination of 30-60 % (Table 1).The soil was classified as silty/sandy loam with depths >40 cm and a humus layer of around 3 cm.In PAM1 and PAM2, canopy closure was 5 to 10 % higher than in SMO1 and SMO2 (82-88 %).The stem density (Fig. S3) ranged between 304 N ha − 1 (PAM2) and 472 N ha − 1 (SMO1).The amount of deadwood varied largely between stands, with lowest values in PAM1 (1.04 m 3 /ha) and highest in SMO2 (28.7 Fig. 1.Location of the four study stands (grey circles) in the Western Rhodopes, southern Bulgaria: stands with "single-tree" selection cuttings are close to Pamporovo; stands with "group-tree" selection cuttings are close to Smolyan.The natural distribution range is shown for spruce (orange), beech (green), and fir (blue) based on distribution maps by Caudullo et al. (2017).Location of the climate station Chepelare is indicated by a black cross and is located at 1500 m a.s.l. in a distance of 5-10 km to the study area.The Walter-Lieth diagram (Walter and Lieth, 1960) characterizes the local climate at the study sites for the period 1990-2019.Map was created using raster data from www.naturalearthdata.com(public domain).m 3 /ha), the latter caused by some wind-thrown trees.Tree height (h) was assessed in dominant and co-dominant fir and spruce trees and ranged between 19 m and 42 m, with spruce being slightly higher on average than fir.The diameter at breast height (dbh) of these trees ranged between 27 cm and 86 cm.A h/dbh-ratio of around 60 and a relative living crown length >50 % indicated high stability and vitality of the studied forest stands.All methods for assessment of soil and stand characteristics as well as the stand descriptions are shown in detail in Schmied (2020).

Management history
According to the forest management plan, the four study stands originated from natural regeneration and were managed traditionally as coniferous-dominated forests mainly for recreational purposes.Therefore, management disadvantaged beech trees and saplings in 1950-2000.However, since the beginning of the 21st century and in the context of climate change, alternative silvicultural methods were applied to enhance the structural heterogeneity of the standsthe growing stock was reduced, and the natural regeneration promoted.The two stands close to Pamporovo (PAM1 and PAM2) were treated by moderate "single-tree" selection cuttings with an intensity of 10 % of the standing biomass, and the two stands close to Smolyan (SMO1 and SMO2) -by more intense "group-tree" selection cuttings with an intensity of 25 % of the standing biomass.These intervention schemes created gaps with a size of 150-250 m 2 in PAM1 and PAM2, and 300-500 m 2 in SMO1 and SMO2, as was estimated during field work.At the sites near Smolyan, wind-throw events have additionally enlarged the canopy gaps; particularly in the centre of SMO2, one large gap of approximately 800 m 2 with a high amount of deadwood was found (Table 1).

Sampling design
In each stand, we established 25 sample plots ("plots") in a grid of 5×5 with a distance of 30 m to each other.The minimum distance to the forest edge was 15 m.This resulted in a rectangular study site of 150 m side length and an overall size of 2.25 ha with 25 systematically organized plots respectively (Fig. S4).Only one site, PAM1, was slightly larger (2.7 ha) and consisted of 30 plots in a 6×5 grid.
A total of 105 plots were sampled in September 2019.Each plot consisted of two concentric circles differing in size: a small plot of 25 m 2 (radius = 2.82 m) and a large plot of 100 m 2 (radius = 5.64 m).

Small-plot measurements
In each small plot, the species-specific regeneration density (RD) and regeneration height were assessed.Regeneration was defined as all individuals ≥10 cm height and <8 cm dbh.Their height was measured with a precision of 1 cm as the vertical distance between the ground surface and the top of the terminal shoot.The number of regeneration plants of spruce, fir, and beech was recorded in three height classes (HC) -HC1: 10 -29 cm height ("seedlings"), HC2: 30 -129 cm height ("saplings"); HC3: 130 cm height − 7.9 cm dbh ("secured regeneration").The regeneration density was calculated for each plot, tree species, and HC as the number of plants per hectare (N ha − 1 ).
Ground cover (%) for the entire small plot area was initially assessed in the field in 16 categories (Table 2).These categories were first distinguished into three surface types ("vegetation", "substrate", and "soil").The surface type "vegetation" comprised nine categories, focusing on vascular plants; the surface type "substrate" subsumed five categories, including mosses, stones, or deadwood (stumps, coarse woody debris, or branches); and the surface type "soil" comprised two categories, litter and mineral soil.In a second step, the categories within each surface type were grouped according to the patterns of variation resulting from a principal component analysis PCA (Fig. S5).This analysis resulted in seven ground cover groups (Table 2), which were verified based on their known ecological effects on regeneration density from the literature (e.g.Balandier et al., 2006;Dyderski et al., 2018;Kalt et al., 2021).These seven ground cover groups were finally used for the statistical analyses.
The browsing status of the apical shoot (damaged/undamaged) was recorded in all individual plants without distinguishing between summer and winter browsing.Plot-specific browsing intensity (BI, %) was then calculated for each tree species and HC as the percent of individuals with browsed apical shoots in the total number of trees in the respective category.
Additionally, in the two tallest and not damaged (i.e., diseased or browsed) individuals per species and plot (N = 446), the last three annual height increments of the apical shoot were measured with a precision of 1 mm.By selecting the tallest plants, we aimed to minimise the effects of competition from the understory and thus get a better reflection of the influence of the overstory (Dȃnescu et al., 2018).This resembles the gap-filler approach (Lertzman, 1992), which is used to assess forest dynamics with a focus on future species composition (e.g., Nopp-Mayr et al., 2020).The averaged individual increments over the last three years (I av ) were then used for statistical analyses.

Large-plot measurements
At the larger plots, we recorded the diameter of all spruce, fir, and beech trees with dbh ≥ 8 cm and calculated the total basal area (BA, m 2 ha − 1 ).Further, the relative BA of fir, spruce, and beech trees (BA fir , BA sp , BA be , respectively) was calculated as a proportion (%) to the total BA of the corresponding plot.This variable was used as a proxy for the presence of mother trees near the small plots.
To analyse the influence of forest canopy openings on RD, the canopy gaps ≥ 100 m 2 were monitored within each study site.Large plots, which overlapped with such canopy gaps, were categorised as "with gap" and "no gap".This variable is further called "gap overlap", and was included in the statistical analysis.
In addition, we determined the type of humus layer (mull/mor) at each plot by excavating an approximately 15 cm thick block of the forest floor.Light conditions (direct site factor DSF, indirect site factor ISF and total site factor TSF) were assessed using a solariscope (Behling, 2015) at two meters height above the plot centre (Annighöfer et al., 2019).The best-fitting image regarding light transmission and canopy cover was selected visually in the field and later controlled in the lab.

Statistical analysis and modelling
Generalised linear mixed models (GLMMs) with a negative binomial distribution were used to analyse the factors driving RD, as they are suitable for non-negative count data with zero-inflation (Zuur et al., 2009).The models were built with the glmmTMB-package (Bolker, 2020).To understand the role of light conditions on I av , general linear regression models (Gamma family with a log link) were applied for each tree species and HC.The effect of key variables on the respective dependent variable was demonstrated using the packages effects (Fox and Weisberg, 2019) and ggeffects (Lüdecke, 2018).All statistical analyses were provided by the software R, version 4.0.0 (R Core Team, 2020).Data exploration was based on the suggestions of Zuur et al. (2010).

Regeneration density modelling
The statistical models used RD per tree species and height class as dependent variables (Table 4).As explanatory variables, measures related to: (1) light, (2) microsite and browsing, (3) management, and (4) ground cover were included.The summary of all explanatory variables is shown in Table 3.As browsing was present only in fir and beech, the explanatory factor "species-specific browsing intensity" was included only in the models of these two species.
As the main interest was on the species-specific effect of DSF, ISF or TSF on RD in different HC, four different versions of the model were formulated: (1) a basic model without including light, (2-4) a model with DSF, ISF (Table 5) or TSF (Table S1) as a light variable.The four model versions were then compared in their performance using AIC (Akaike information criterion; Akaike, 1974) with the function AICtab (Bolker and R Development Core Team, 2020), and the model with the lowest AIC is presented in the result section.The goodness-of-fit of all models was checked with DHARMa package in R (Hartig, 2020) for dispersion, zero-inflation, and outliers.Further model diagnostics were applied with the package performance (Lüdecke et al., 2021) for collinearity, zero-inflation, and the calculation of R 2 .Pairwise comparisons between groups (Wilcoxon rank-sum tests) were applied for differences in the effect of "gap overlap" and "treatment" on RD (Table S2).

Height increment modelling
Average height growth, I av , served as the dependent variable.The explanatory variables were "species" and the "light factor" (DSF, ISF or TSF).Interactions between "species" and "gap overlap" were investigated by pairwise comparisons with Wilcoxon rank-sum tests and Bonferroni correction between groups.For each HC, three model versions differing in the "light factor" were formulated: the results differed only marginal in the effect size, so we used the TSF as it accounts for both, direct and diffuse light, and therefore best represents PPFD.The distribution of I av data was previously tested by the descdist function of the fitdistrplus package in R (Delignette-Muller and Dutang, 2015).Outliers with Cook's distance > 6/N (N = dimension of the dataset) were excluded from the analysis, and the models were then refitted (Crawley, 2009).Outliers' exclusion (n ≤ 3 per model) did not change the model outcome but did slightly improve the model fit.Fitted models of beta (betareg package in R; Cribari-Neto and Zeileis, 2010; Grün et al., 2012) and gamma distribution were compared via posterior predictive checking (Gelman et al., 2021;Lüdecke et al., 2021) and residual diagnostics (Hartig, 2020).The models with gamma distribution showed a better fit.

Stand characteristics
The total BA ranged between 36 m 2 /ha (in PAM2) and 54 m 2 /ha (in SMO2).Fir had the highest proportion in the overstory (44 % and 50 %), apart from SMO2, where spruce dominated (49 %; Table 3).The proportion of beech was higher in PAM1 and PAM2 (25 % and 35 %) than at the Smoljan sites (7 % and 15 %), but generally had a lower proportion than fir or spruce, except for one site, where spruce had the lowest BA (SMO2, 14 %).
The number of plots overlapping with gaps ranged between 7 and 21
Site Species Regeneration density (trees/ha) D. Ambs et al.
per site with a higher number and larger gaps at both Smoljan sites (n = 16 and 21) than at the Pamporovo sites with smaller-sized gaps (n = 7 and 13; Table 3).

Light and microsite characteristics
The light conditions on the sites PAM1 and PAM2 were very similar, with a median of 2 % for DSF and 3 % for ISF (Fig. 2).In contrast, SMO1 and SMO2 had a median of 4 % and 9 % for DSF, and 9 % and 12 % for ISF (Fig. 2).Among all plots, the greatest variability in light conditions was found at the "group-tree" treatment sites near Smolyan for both ISF and DSF (0 % to 43 %; Fig. 2).
Across all sites and height classes, 5 % to 10 % of the fir trees were browsed, while beech showed the highest proportion of browsed trees (15 % to 30 %; Table 3).The strongest impact was in both species on HC 2 (fir: 13%; beech: 45%), whereas in HC1 and HC3 browsing intensity was 6 and 4% (in fir) and 10 and 16% (in beech), respectively.In spruce, browsing was negligible (<1%), and was therefore not included in further analyses.

Regeneration density
Regeneration plants were found on all 105 plots: fir was present on 98 % of the plots, followed by spruce (61 %) and beech (48 %).The maximum RD was found for all species in HC1, with fir reaching 55000 N ha − 1 , followed by beech (22500 N ha − 1 ), and spruce (14000 N ha − 1 ).
In general, the RD of seedlings was highest and decreased with increasing size of the regeneration plants (Table 4).Across all sites, fir had the highest mean RD in HC1 (7000 to 14000 N ha − 1 ) and in HC2 (4000 to 9000 N ha − 1 ) compared to the respective HC of spruce and beech.However, in the secured regeneration class, fir was outcompeted by beech in PAM 1 and PAM2 (120 to 300 N ha − 1 vs. 400 to 600 N ha − 1 ; Table 4).Silvicultural approaches of "single-tree" cutting were pursued

Table 5
Summaries of the GLMMs per tree species (fir, spruce, beech) and height class HC (HC1 = 10 -29 cm, HC2 = 30 -129 cm, HC3 = 130 cm height -7.9 cm dbh) with coefficient means ± standard error.Model type describes the best fitting version of the four models tested.Distribution shows the specification of the negative binomial distribution, and R 2 represents the squared correlation between the model's actual and predicted response.at these sites.Spruce had the lowest RD in all height classes at these sites, particularly in the secured regeneration class HC3, where spruce was rare (40 N ha − 1 in PAM1) or absent (as in PAM2).The "group-tree" selection approach (at SMO1 and SMO2) led to average spruce RD in all HC.At these sites, beech regeneration plants were rarely found in HC1 and HC2 (0 to 200 N ha − 1 ) but had relatively high density in the secured regeneration class (150 to 350 N ha − 1 ).

Fir
In HC1, the basic model without including light showed the best performance and explained 33 % of the variability (R 2 = 0.33; Table 5).From the microsite-related variables, none of the predictors were suggested to be detrimental for RD of the fir seedlings.Treatment also showed no impact on RD of fir seedlings (Fig. 3).From the "ground cover" variables, TallBerry, MossStone, and Litter had a significantly negative, and ShrubsGrass a significantly positive effect.The most substantial negative effect on HC1 of fir was detected on positions covered by TallBerry, where at 40 % cover degree, the RD of fir saplings was halved (Fig. 4).In HC2, the best model performance was noticed when including ISF (R 2 = 0.60; Table 5).A higher RD of fir was associated with a higher proportion of ISF, while the proximity to gaps ("gap overlap") had a significantly negative and strong effect (-0.39; p < 0.05).From the "ground cover" variables, a higher amount of TallBerry (Fig. 4) and ShrubsGrass affected RD saplings negatively.The substrate types of Litter and CWD were also negatively influencing RD of HC2 in fir.The best-performing model for HC3 included ISF; the "proximity to gaps" was associated with a strong significant increase in RD, while a high BA was negatively influencing RD.

Spruce
The basic model fitted best for the spruce seedlings (HC1), but due to the low explanatory power of the model (R 2 = 0.09; Table 5), we will not interpret it further.The best model for HC2 contained DSF and had the highest explanatory power (R 2 = 0.40) among all other models.A high BA sp was detrimental for RD, but the effect size was small (− 0.01; Table 5).The "proximity to gaps" had the strongest positive effect (0.85; Table 5).A high cover degree of ShrubsGrass, Herbs and CWD was associated negatively with RD; however, these variables had a small effect and were insignificant (p ≤ 0.1; Table 5).The best-performing model for HC3 revealed a positive effect of the "group-tree" selection treatment on RD (Table 5; Fig. 3); this effect had a large effect size but was not significant.BA sp was also negatively influencing RD, but this effect was only small (− 0.01; Table 5).In addition, less secured regeneration plants of spruce were found on positions covered with Moss-Stone.Positive relationships were found between RD and the ground cover categories TallBerry, woody debris (CWD), and Litter.

Beech
The best-performing model in HC1 revealed that more direct light increased the density of beech seedlings.The management variables "BA be" and "proximity to gaps" were related to high RD, whereas treatment ("group-tree") had a large negative effect (− 7.54; Table 5; Fig. 3).From the "ground cover" variables, TallBerry was negatively related to RD of beech seedlings, whereas Litter was related positively; however, the effect sizes of both variables were small (Table 5).
The basic model was the best in HC2."Group-tree" selection treatment had a large negative effect on RD in HC2 (− 3.47; Table 5; Fig. 3).Lower RD of beech saplings were associated with moder and presence of TallBerry vegetation type.
The basic model for HC3 showed that higher RD was related to a smaller proximity to gaps (0.88; Table 5) and moder as a substrate type.However, the explanatory power of HC3-model was very low (R 2 = 0.06).

Descriptive analyses
Height increment I av showed only little differences between the tree species in HC1 and was around 1.40 cm (Table 6).The increment of HC2 saplings was highest in beech (ca.6.5 cm), followed by fir (ca. 5 cm), and spruce (ca. 4 cm).Across all three tree species, the annual height increment of secured regeneration (HC3) was three to four times higher than those of HC2 and nearly ten times higher than the increment of HC1 seedlings.While fir and beech (HC3) increased in height around 14.0 cm per year, the annual increase in spruce was about 9.5 cm (Table 6 and Table S3).
The proximity to a gap was important for the saplings (HC2) where all species showed highest I av when growing on locations close to gaps (Fig. 5).Beech saplings seemed to profit stronger from the gaps as compared to fir and spruce (Fig. 5) and showed a twice higher I av compared to saplings growing on locations without gaps.In HC3, due to the low number of trees under closed canopy the test results will not interpreted here (Table S2).Fig. 3. Treatment effects ("single", "group") on regeneration density separated by height classes (HC1 = 10 -29 cm "seedlings", HC2 = 30 -129 cm "saplings", HC3 = 130 height -7.9 cm dbh "secured regeneration") for fir (blue), spruce (orange), and beech (green).Black dots with error bars indicate the predictions of the effects; the raw data are shown in the species-specific colour.All continuous variables were held constant at their mean; all categorical variables were weighted proportional to their sample size (Fox and Weisberg, 2018).Note scaling differences on the y-axis.Fig. 4. Effects of the cover degree of TallBerry on regeneration density (RD) of fir separated by height class (HC1 = 10 -29 cm "seedlings", HC2 = 30 -129 cm "saplings", HC3 = 130 -7.9 cm dbh "secured regeneration").The 95 % confidence intervals are shown in the species-specific colour.All continuous variables were held constant at their mean; all categorical variables were weighted proportional to their sample size (Fox and Weisberg, 2018).Note scaling differences on the y-axis.

Factors driving height increment
Since the I av models led to similar results when DSF, ISF or TSF were included, we only show the results of modelling with the total site factor (TSF).A general positive effect of TSF on I av was found in HC2 and HC3, whereas I av of the seedlings (HC1) was light-independent (Table 7; Fig. 6).
We found species-specific differences, but not in all height classes.For the seedlings (HC1), fir showed the highest annual increment (1.48 cm; Table 7), followed by spruce (1.23 cm).In contrast, beech seedlings showed a lower I av (1.00 cm; Table 7).In HC2, beech had the largest increments (5.38 cm), followed by fir (3.56 cm) and spruce (2.79 cm).The height increment of the secured regeneration was similar in fir and beech (8.90 cm and 9.89 cm, correspondingly) and significantly lower in spruce (5.70 cm).

Discussion
In this study, we investigated the main factors controlling regeneration dynamics in four mixed stands from the Western Rhodopes undergoing a conversion from formerly spruce-dominated forests to increasingly mixed stands.Our study (1) focused on the main drivers influencing tree regeneration in the initial process of conversion and (2) evaluated two management scenarios with varying harvesting intensity.

Fir
In all study stands, fir is the species with the highest regeneration density in the seedling and sapling stage of development.However, in the secured regeneration class, fir had the highest RD only at the two SMO sites, which were characterized by large gaps and more available light.Fir has been considered to be established well under relatively closed canopies, whereas under direct-light conditions the recruitment often was impeded due to its lower competitiveness compared to competing ground vegetation (Grassi and Bagnaresi, 2001;Klopčič et al., 2015;Mosandl and El Kateb, 1988;Vencurik et al., 2020).Our results are in line with these studies, as fir saplings had lower RD at plots located in proximity to gaps.In addition, we detected a negative effect of TallBerry on fir seedling and sapling density, most likely due to competition (see also Thrippleton et al., 2016).Fir saplings are known to be relatively sensitive to late frosts (Kohnle et al., 2011) that occur more often in gaps or open areas than under closed canopy.The positive effect of indirect light on fir saplings' RD that we found in the present study reflects typical low PPFD conditions in stands with small to medium canopy openings and is considered to be beneficial for fir seedlings' ingrowth into the sapling stage (Mosandl and El Kateb, 1988;Orman et al., 2021).In line with this, our results showed that fir seedlings and saplings are susceptible to be located in canopy openings; in contrast, secured regeneration had higher densities at these sites.However, we did not find a treatment effect, which indicates that the change in key microsite factors close to gaps (e.g., light supply, cover degree by Tall-Berry vegetation type) plays a stronger role in RD of fir seedlings and saplings than the factors related to the management type itself (e.g., management-induced mortality/disturbance of plants during logging).The secured regeneration of fir was mostly driven by gap proximity, which may include the influence of a management-related artefact, i.e., the forester's decision to create and/or enlarge existing gaps close to or around hotspots with already existing, mostly secured regeneration.Fig. 5. Height increments I av (cm) of individual regeneration plants for the period 2017-2019 in plots "no gap" (grey) and "with gap" (coloured: fir = blue, spruce = orange, beech = green).Boxplots sharing a letter show no significant difference in I av between increment of plants on "no gap" and "with gap" plots (Wilcoxon tests, p < 0.05).Please consider the different scale on the y-axis.

Spruce
In our study, the RD of spruce remained relatively low in all height classes.Norway spruce is a tree species known for its high morphological and physiological plasticity in reaction to increasing light (Grassi and Bagnaresi, 2001).The relatively frequent seed production, the long distances of seed dispersal, as well as the low susceptibility to late frosts help spruce to colonize canopy gaps faster as beech and fir (Gratzer and Waagepetersen, 2018;Paluch et al., 2019).Additionally, spruce has been shown to benefit from direct radiation (Bednář et al., 2022), while closed canopies are likely to impede seedling establishment (Mosandl and El Kateb, 1988).Although the statistical models showed not always a significant influence of proximity to gaps (and/or increasing light) on RD of spruce, the inclusion of light as predictor variable improved model performances for HC2 and HC3, indicating that light-related microsite factors (e.g., temperature, water-balance) might be an important factor for spruce establishment on these mountain forest sites.The proximity to gaps was positively related for spruce saplings RD, but not for secured regeneration.This outcome may be explained due to the low abundance of secured regeneration plants, which made our spruce HC3-models noninterpretable.Notable is the significant negative influence of increasing proportion of spruce trees in the overstory which is contrary to previous findings of Lundqvist and Fridman (1996) from Boreal forests.In a study in the Eastern Alps, Diaci et al. (2020) point out a possible below-ground competition for water between recruitment and mature trees, particularly during prolonged summer droughts.This is supported by recent findings by Simon et al. (2019) for conifers but not for beech, which may be attributed to higher stem flow in beech trees during precipitation events (Patzner, 2004;Zirlewagen and von Wilpert, 2001).During extreme years with warm and dry late summers, which are characteristic at our study sites, a high proportion of mature spruce trees may have led to high mortality of shallow-rooting spruce regeneration plants, thus, creating competitive advantages for fir and beech recruitment.In colder habitats, Brang (1998) observed that spruce seedlings on north-oriented slopes preferentially grow on warmer microsites with more direct light.However, competing vegetation (i.e., TallBerry) is becoming a problem under enhanced direct light.To overcome this competition, spruce recruitment often occurs on elevated microsites or on coarse decaying wood debris (Dountchev and Zhelev, 2015;Hunziker and Brang, 2005;Vencurik et al., 2020).Our study did not support the importance of woody debris for spruce establishment, which may be due to the low presence of dead wood in an appropriate stage of decay.Most stumps and nurse logs were estimated to be approximately 10-15 years old.As reported by Zielonka (2006) and Tsvetanov et al. (2018), several decades of decay are required until dead wood becomes suitable for the colonization by spruce.Additional negative effects may have been provoked by occasional late summer droughts, causing the nurse logs to dry out and make them less suitable for seedling growth.Therefore, a lack of suitable microsites for early spruce establishment may have restricted RD of spruce seedlings in the study sites in Western Rhodopes.

Beech
European beech shows the best regeneration success under moderate shade conditions, which may result from silvicultural operations that create gaps of 100-500 m 2 in size or after small-to medium-scale natural Fig. 6.Visualization of the effects of the averaged height increment (I av ) models for the variable TSF (total site factor) separated by height class HC (HC1 = 10 -29 cm "seedlings", HC2 = 30 -129 cm "saplings", HC3 = 130 -7.9 cm dbh "secured regeneration") and tree species: fir (blue), spruce (orange) and beech (green).The 95 % confidence intervals are shown in the species-specific colour.Please consider the different scale on the y-axis.
D. Ambs et al.
disturbances (Feldmann et al., 2018;Hobi et al., 2015).In gaps larger than 500 m 2 , the competitiveness of beech seedlings may decrease (Petrovska et al., 2023), resulting in a less successful establishment (Bílek et al., 2013;Kelemen et al., 2012;Pretzsch et al., 2015).We observed a relatively high density of beech seedlings at sites with a relatively closed canopy following "single-tree" cutting approaches (PMO1, PMO2), which also favour shade-tolerant species like fir (Käber et al., 2021).This is in line with the low mortality rates reported for beech recruitment under deep shade conditions (Petrovska et al., 2022).Increased presence of TallBerry vegetation had a negative effect on beech regeneration density, indicating that opening the canopy too much might be detrimental in the long-term.
The presence of mature beech trees increased regeneration density, which can be explained by the limited distance of seed dispersal in beech (Gratzer and Waagepetersen, 2018;Paluch et al., 2019).This effect indicates that beech recruitment may have been initially limited by the relatively low presence of mother trees in the studied stands, particularly in SMO1 and SMO2.In the seedling stage of development, a higher stem density of beech compared to fir has been observed at microsites with high light availability, whereas higher densities of fir were observed under relatively closed canopy (see Bottero et al., 2011).In the late sapling and secured regeneration stage, Bottero et al. (2011) reported higher shares of beech under more shady conditions, which is in line with our observations from the sites with "single-tree" selective cuttings.

Overall browsing intensity
Damage caused by game animals can be a challenge for regeneration processes in mixed mountain forests and is a widespread disturbance factor in fir-spruce-beech stands, especially in the Alpine region (Liss, 1988;Prietzel and Ammer, 2008;Simon et al., 2019).In the Western Rhodopes, due to the high presence of wolfs and bears (Genov and Kostava, 1993), we expected less browsing intensity on the regeneration, similar to other forests in south-eastern Europe (Bottero et al., 2011;Tsvetanov et al., 2018).In our study, we found stronger browsing preference for fir and beech than for spruce, which is in accordance to reports for mixed forests in the Alps (Kupferschmid et al., 2019;Liss, 1988).For fir the browsing intensity in HC1 and HC3 was relatively negligible.The damage on HC2 with 13 % was similarly high as reported by Kupferschmid (2018) for beech-fir forests in Switzerland, but much lower than the damage recently reported by Thom et al. (2022) (around 50 to 60 %) for forest sites in a German national park.However, due to the very high stem density of the regeneration, overall damage is less critical to the successful regeneration of these stands (e.g., Nopp-Mayr et al., 2020).To our knowledge, this is the first study reporting on browsing damage in Bulgarian forests; therefore, further examinations are required to analyse the effect of ungulates on species composition in these forest regions.

Main factors driving height increment
Canopy openings and light are known as crucial drivers of forest regeneration growth dynamics (e.g., Annighöfer, 2018;Dȃnescu et al., 2018;Diaci et al., 2012;Preuhsler, 1981;Stancioiu and O'Hara, 2006;Thom et al., 2022).In the present study, the effect of light was more pronounced with increasing tree size, which is confirmed by findings of Schmid et al. (2021).Additionally, we detected that differences related to tree ontogeny were species-specific: in HC1, fir increment was most sensitive to TSF, followed by spruce and beech; in HC2 and HC3, beech was most light-sensitive, followed by fir and spruce.Other studies also demonstrated that the shade tolerance of a species can change during tree ontogeny (e.g., Annighöfer et al., 2017;Käber et al., 2021;Niinemets, 2006).In all height classes, however, spruce showed the lowest I av compared to beech and fir, high-lighting spruce as the least shadetolerant tree species among all tested (see also Käber et al., 2021).
In the seedling stage, no significant dependence of I av to increasing light conditions could be detected, in the seedling stage.This suggests that other factors, such as micro-relief or topsoil moisture, limit growth more than light at this early developmental phase (Bebre et al., 2020;de Frenne et al., 2021;McCarthy, 2001).However, height increment under canopy showed no significant difference to those in proximity to gaps.One reason for these ambiguous tree growth responses could be that we did not further differentiate gap shape or size (Diaci et al., 2020;Simon et al., 2019).Another explanation could be related to additional resource allocation to root biomass for better water access at sites experiencing drought episodes (Nikolova et al., 2020).Saplings and secured regeneration of fir and beech responded more strongly in their height increment to increasing radiation than spruce.One reason could be the insufficient light intensity in the stands of our study (4-15 % TSF).
According to Stancioiu and O'Hara (2006), the increment of spruce is smaller compared with fir and beech under low light conditions (<35 % of open canopy).This competitive disadvantage may lead to lower survival rates (Kobe et al., 1995;Löf et al., 2007) and reflects the lower shade-tolerance of spruce compared to fir and beech (Diaci et al., 2020;Käber et al., 2021;Vencurik et al., 2020).These findings are also supported by Laiho et al. (2014), who emphasize the importance of the enlargement of canopy openings for taller spruce saplings.In a different direction, Unkule et al. (2022) argued that high summer temperatures had a negative effect on spruce and beech saplings, while fir was not affected.Therefore, the present late-summer droughts (see Bocheva et al., 2009), which can be exacerbated by direct light in proximity to canopy openings, could explain for the relatively low height growth of spruce.Hence, the low densities of spruce may result from its relatively poor performance at these marginal site conditions, which increases the risk of spruce to become less present or even disappear in these forests.
In contrast, fir and beech I av increased sharply with light in the larger size classes, indicating, in both tree species, comparable vigour and optimal growing conditions for regeneration development.One factor which may influence the species growth reaction can be the speciesspecific acclimatization ability of fir and beech after rapid exposure to increased radiation ( Čater and Diaci, 2017).The limited ability of fir saplings to exploit high-light conditions may be a competitive disadvantage in large canopy gaps and thus may restrict the recruitment processes of fir to relatively small gaps (Grassi and Bagnaresi, 2001).

Conclusions and management recommendations
Our case study was based on data from four forest sites that differed in management and revealed an overall high recruitment potential in the Western Rhodopes.However, studied natural regeneration processes indicated a possible change in future forest composition with spruce becoming less present in the studied stands.The results showed that the "single-tree" selection approach tended to promote fir dominance, while the "group-tree" selection approach promoted beech and fir, but also spruce, especially when suitable microhabitats were available.The growing conditions created by larger gaps (300-1000 m 2 ) seemed to successfully convert the formerly spruce-fir dominated forests to climate-adapted fir-beech-(spruce)-mixed stands, which is in line with previous studies (Brang et al., 2014;Lafond et al., 2017).In contrast, when creating smaller gaps, beech and fir will dominate future forest composition.If ungulate density remains under control, browsing did not seem to be a crucial disturbance factor in the forest dynamics of the Western Rhodopes.Since rejuvenating forest stands in a "close-to-nature approach" needs several decades, foresters may not act schematically or one-dimensional, but case by case in a staggered/patchy way.Moreover, management approaches should focus on constraining competing ground vegetation, like ferns, raspberries, blackberries, and tall forbs, by slowly and steadily opening the canopy to increase structural heterogeneity.This should ensure a successful ingrowth of differing size classes of regeneration, and thus the conversion to mixed fir-beech-spruce stands.Appropriate silvicultural approaches should be transferred into prescriptions for forest practitioners (Pretzsch and Zenner, 2017).

D. Ambs et al.
However, further research on gap dynamics and microclimate under different management and/or disturbance regimes is required (Thom et al., 2022).
In conclusion, we revealed that less intense and more frequent, punctual interventions following advanced regeneration (irregular shelterwood or selection systems) create more optimal conditions for tree recruitment than regular harvests with greater intensity (regular shelterwood).The former also creates more resilient forest structures under conditions of intense natural disturbances (Dountchev and Zhelev, 2016) and better habitat quality for many species (Nikolov et al., 2022).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Boxplots of direct site factor (DSF) and indirect site factor (ISF) by study site (PAM1, PAM2, SMO1 and SMO2).Different letters refer to pairwise comparisons of DSF (lowercase and ISF (capital letters) among the sites (non-parametric Wilcoxon rank sum tests with Bonferroni correction).Means not sharing the same letter are significantly different.TSF data are not shown.

Table 1
Soil and stand characteristics of the stands PAM1, PAM2, SMO1 and SMO2.H top5 reflects the averaged tree height of five tallest trees per tree species (fir, spruce) and site.Data shown as means ± standard error.

Table 2
List of the 16 initial ground cover variables which served for the construction of the seven ground cover groups used for the statistical analyses.

Table 3
Summary statistics of the input variables considered in the statistical analyses of regeneration density.Continuous variables are shown by means ± standard error.The categories set as a reference in the statistical models are shown in bold italics.
D. Ambs et al.