E ﬀ ect of Invasive Rhododendron ponticum L. on Natural Regeneration and Structure of Fagus orientalis Lipsky Forests in the Black Sea Region

: Biological invasions threaten global biodiversity and forest ecosystems; therefore, it is necessary to use appropriate strategies for combating the spread of invasive species. Natural regeneration of eastern beech ( Fagus orientalis Lipsky) is considerably limited by an aggressive invasive shrub, pontic rhododendron ( Rhododendron ponticum L.), in the Black Sea Region of Turkey. Therefore, the future character of the region’s forests is uncertain. The aim of this research was to evaluate the structure of beech forests with di ﬀ erent management regimes of rhododedron and to determine the interaction among tree layer, rhododendron cover, and natural regeneration in Düzce Province using the FieldMap technology. The following variants of forests were compared: without intervention (control) and three and six years after rhododendron clearance. The results showed that tree density ranged between 175–381 trees ha − 1 and stand volume between 331–589 m 3 ha − 1 . The horizontal structure of the tree layer was mostly random, and the spatial pattern of natural regeneration was aggregated. Recruit density and height in the beech stands were signiﬁcantly di ﬀ erentiated due to the inﬂuence of presence or absence of invasive rhododendron. Rhododendron cover ranged between 81%–97%, and woody stems amounted to 72,178–86,884 ha − 1 in unmanaged forests. Canopy in the overstory did not have a signiﬁcant e ﬀ ect on the density of regeneration and rhododendron cover. Tree layer had a signiﬁcant negative inﬂuence on natural regeneration within a 4 m radius on the plots without rhododendron. However, on the plots with dense rhododendron cover, tree layer had a positive inﬂuence on regeneration within a 1.5 m radius. Natural regeneration density was signiﬁcantly higher when rhododendron was cleared than the plots without intervention. On the plots without woody clearance, there was an insu ﬃ cient regeneration (113–619 recruits ha − 1 ); however, they had higher mean height compared to the sites without rhododendron. After three and six years of rhododendron clearance, the numbers of recruits in natural regeneration were 63,981 ha − 1 and 105,075 ha − 1 , respectively. In conclusion, invasive spread of rhododendron was a limiting factor of the prosperous regeneration and tree species diversity, and manual clearance of rhododendron is recommended in managed beech forests of the study region.


Introduction
In the past, management of forests in Turkey was aimed at maximizing timber production regardless of sustainability of the multiple forest functions. However, forest stands did not receive any silvicultural treatment, all standing trees are above their economic rotation age, and almost half of the Turkish forests are considered to be degraded in terms of wood production [1]. The majority of Decisions on the methods of differentiated management of forests on ecological principles are complicated to ensure their ecological stability and biodiversity [28,56]. Models can be applied for support of a decision-making and planning process that helps finding the optimum solution to forest management [57]. Forest management and silviculture treatments (including thinning, promoting optimal growing conditions, reducing hosts, soil preparation, and consistent monitoring) then play the main role to control invasive alien species [58,59] and protect, in our case, threatened beech forests from aggressive invasive pontic rhododendron [60]. However, optimal management strategies should always be based on in-depth investigation and evaluation of species-specific site and stand conditions, especially stand structure and development of the natural and near-natural forests [61]. Our study can contribute to the solution of this serious problem.
A general objective of our study was to evaluate effect of invasive pontic rhododendron on the structure of eastern beech stands in the north-west of Turkey in the Black Sea Region, with an emphasis on natural regeneration. Specific objectives were (1) to evaluate the parameters describing production, growth, diversity and structures of eastern beech stands, (2) to quantify and evaluate the effects of invasive pontic rhododendron and tree layer on the prosperity of natural regeneration in commercial and noncommercial forests, (3) to assess the effects of stand structure on the occurrence and expansion of pontic rhododendron, and (4) to quantify and evaluate the interactions among rhododendron, natural regeneration, stand production, and structural diversity.

Study Area
This study was conducted in the close-to-natural forest stands with dominant eastern beech in the Province of Düzce in the Black Sea Region (BSR) of Turkey. Five permanent research plots (PRP) were established in the mountain regions 14 km northeast of Düzce in 2016. Three PRP (1, 2, and 5) with historical minimal harvesting operations were situated in the stands with different tree canopy closure in the overstory and dense rhododendron presence in the understory. On PRP 3 and 4, rhododendron was removed manually by grabbing the roots three and six years ago, respectively, and left to decompose in soil. Moderate thinning of the tree layer was carried out on these two PRPs during rhododendron clearance. The proportion of stand volume logging were 8% (35 m 3 ha −1 ) on PRP 3 and 12% (65 m 3 ha −1 ) on PRP 4. The PRP were located in similar sites in close proximity to each other so that they were comparable ( Figure 1, Table 1). The study area falls in the jurisdiction of the Düzce Forest Management Directorate of Bolu Regional Forest Directorate in the western BSR, which is in the Euro-Siberian Floristic Region [62]. The territory of study area belongs to temperate oceanic climate characterized by a warm dry summer and cool dry winter with a relatively narrow annual temperature range [5]. The mean annual temperature in the area is 13.4 • C and mean temperature of the coldest and warmest months reaches 3.8 • C (January) and 22.6 • C (July), respectively. Amount of annual precipitation ranges around 820 mm with maximum in December (100 mm), mainly in the form of snow. The number of rainy days is 135 (Turkish State Meteorological Service). Soil texture ranges from clay and clay loam to loamy clay. The forest floor surface is 3-5 cm deep. The A-horizon is 10 cm deep, underlain by B-horizon of 40-60 cm in thickness. Stoniness is medium, ranging from 10% to 30% by volume. Soil profile depth is 110 cm, though rooting depth may go deeper between rocks. Soils in this region of Turkey have been placed in the Inceptisols order, under the Umbrepts suborder, and are classified as Typic Haplumbrepts according to the USDA soil taxonomy [62].

Data Collection
The Field-Map technology [63], which is a hardware tool, was used for establishment of five PRP of 40 × 40 m in size (0.16 ha). This hardware tool helped to determine the tree layer structure on each PRP and located the position of all individuals with breast height diameter (DBH; at 1.3 m above ground) >4 cm. The crown projection area of each individual, minimally at four directions perpendicular to each other, were determined. The crown measurements were used to determine the spatial pattern and canopy cover. Diameter, height, and height to live crown base were also measured for each tree layer. Diameters of the tree layer were measured with a Blue Mantax metal calliper (Haglöf, Sweden) to the nearest 1 mm, and heights were measured with a Vertex laser hypsometer (Haglöf, Sweden) to the nearest 0.1 m.
On each PRP, the position, height, and crown width were measured for recruits with height ≥100 cm and DBH <4 cm. Natural regeneration (recruit) from one year of age (DBH <4 cm) was measured on four nested subplots of 10 × 10 m in size within each PRP ( Figure 2). The position, height, and crown width were recorded for each recruit. The crown width of recruit was measured using a height-measuring pole to the nearest 1 cm.
Basic dendrometric characteristics of pontic rhododendron estimated on each subplot includes cover percentage, stem number, basal stem diameter, and height of individuals. In the total coverage percentage of the rhododendron, each plot was divided into the 100 grids (4 × 4 m grids) by Field-Map technology where cover percentage was subjectively determined for each grid (accuracy 5%). The number of all stems and basal diameter of stems at 10 cm above ground were measured by a metal caliper (accuracy 1 mm), while height was determined by height-measuring pole (accuracy 5 cm) on two subplots, which were selected randomly using the random numbers generated (RNG) in MS Excel.

Data Analysis
Structural and growth parameters, horizontal and vertical structure, production potential of sites, and biodiversity were evaluated on each PRP. The stem volume of eastern beech trees was calculated using the previously developed tree volume equations [64] and that of the hornbeam and oak trees using tree volume equations developed by Petráš and Pajtík [65]. A list of equations used for tree stem volume calculations is given in the Supplementary Files. The curves of height-diameter relationship were constructed using the Näslund function [66]. The fitting performance of the height-diameter model was evaluated using the coefficient of determination (R 2 ). Variances in stand parameters were indicated by standard error of the mean (SE).
In order to determine the spatial distribution pattern of the individuals on each PRP, Pielou-Mountford index [67,68], Clark-Evans index [69], and Ripley's L-function [70] were computed. David-Moore index [71] was used to determine the distribution indices based on tree frequency in particular quadrats. A list of formulae used for determining the spatial distribution patterns is given in the Supplementary Files. The quadrat size of 10 × 10 m was chosen in the overstory and 2 × 2 m in regeneration. The PointPro 2.2 programme (© 2020 Zahradnik) was used for computation of characteristics describing the horizontal distribution of individuals across the research plot. A statistical test on the deviations from the values expected for the random point patterns was performed using the Monte Carlo simulation [72]. For each investigated research plot, 1999 simulations of the Poisson process with the same stand density were performed (α = 0.05, k = 100). The corresponding expected values of these indices were computed by means of numerical simulations for each particular case separately. The confidence limits were estimated using 95% confidence interval.
Evaluation of biodiversity of each research plot involved the computation of a number of indices, such as species richness index [73], species diversity index [74], species evenness index [75], diameter differentiation index and height differentiation index [76], crown differentiation index [77], Arten-profile index [78], and total diversity index [77]. A complete list of formulae used for calculation of biodiversity indices is given in the Supplementary Files. Table 2 shows the criteria of structural indices. The relationship between tree layer and natural regeneration was calculated by the pair correlation function using the R 3.4. software (© 2020 The R Foundation). Stand density index (SDI; [79]-Equation (15)), crown closure (CC; [80]-Equation (16)), and crown projection areas (CPA; sum of measured crown projection per hectare) were also determined for analysis of the horizontal structure on each PRP. Situational maps were created in the ArcGIS programme (© 2020-2010 Esri).
All statistical analyses were carried out using the Statistica 13.5 software (© 2020-2019 StatSoft). Data were log transformed to acquire normal distribution, which was tested using the Kolmogorov-Smirnov test. One-way analysis of variance (ANOVA) and post-hoc Tukey's HSD tests were used to evaluate the research plot-level differences. In addition, relationships among natural regeneration, tree layer, and rhododendron cover were evaluated using the Pearson correlation coefficients. Unconstrained principal component analysis (PCA) in CANOCO for the Windows 4.5 programme [81] was used to analyze the relationships between invasive rhododendron, stand and natural regeneration parameters, and structural diversity. Data were log-transformed, centered, and standardized in this analysis. The results from the PCA were visualized in the form of an ordination diagram.

Growth Parameters and Stand Structure
The number of trees ranged between 175-381 trees ha −1 on PRP 1-5. Eastern beech accounted for 98-100% of the whole stand with the remaining made up by European hornbeam (Carpinus betulus L.) and sessile oak (Quercus petraea Matt.). Stand density index differed between 0.48-0.75 and CPA in the range of 1.90-4.42 ha. Mean stand volume was 427 m 3 ha −1 (± 43 SE) with minimum on PRP 2 (331 m 3 ha −1 ) and maximum on PRP 1 (589 m 3 ha −1 ). Periodic annual increment of the stand ranged between 4.8 and 7.5 m 3 ha −1 y −1 , and mean annual increment varied between 3.3-5.9 m 3 ha −1 y −1 . Average stand basal area was 33.8 m 2 ha −1 (± 3.2 SE; Table 3). Table 3. Structural characteristics of forest stands on permanent research plots in 2016.

PRP
Age Notes: PRP-permanent research plot, age-average stand age, DBH-mean quadratic breast height diameter, h-mean height, v-average tree volume, N-number of trees per hectare, BA-basal area, V-stand volume, PAI-periodic annual increment; MAI-mean annual increment, CC-canopy closure, CPA-crown projection area, SDI-stand density index.
Frequency distribution of diameters on the tree layer and the relationship between DBH and height smoothed by the Näslund function are shown in Figure 3. This function described a large part of the height-diameter relationship on each PRP (R 2 = 0.86-0.94). Diameter distribution shows a typical two-tree layer structure with the main amplitude in DBH class 32-36 cm (max. on PRP 5-63 trees ha −1 ). Based on the height curves, it is significant that the highest trees were located on the PRP 5 lying in a small valley.

Tree Layer Biodiversity
All biodiversity indices and species richness for each research plot are presented in Table 4. Species richness showed low diversity (D = 0.053-0.168). Species heterogeneity and evenness reached minimum values (H = 0.000-0.045, E = 0.000-0.149). The vertical structure according to the A index was medium to strongly diversified (A = 0.414-0.627). Diameter differentiation of the structure was medium (TM d = 0.305-0.452). Height differentiation was mostly low (TM h = 0.203-0.277), only on PRP 5 it was medium (TM h = 0.305). The crown differentiation reached mean values (K = 1.574-1.951). Total diversity for all PRPs indicates an even structure (B = 5.113-5.810). Stands with rhododendron clearance and moderate thinning of tree layer (PRP 3 and 4) had smaller crown differentiation and total diversity compared to other PRP, but there was no difference in other indicators.
The spatial patterns of tree layer, dead wood, and natural regeneration for each PRP are presented in Figure 2. Based on the computed structural indices, horizontal structure of the tree layer was indistinctly aggregated or moderately regular, while random distribution prevailed (Table 4). Research plot 3 was an exception, where the spatial pattern according to R and ICS indices was significantly aggregated (p < 0.05), respectively random according to α index. The L-function showed mostly a random distribution, and it was aggregated within 2 m on PRP 3. It was also aggregated within 1 m on PRP 4, and there was a regular pattern at spacing of 4.5-6 m on PRP 1.

Characteristics of Rhododendron ponticum Understory
Basic dendrometric characteristics of Rhododendron ponticum on PRP 1-5 are shown in Table 5. The rhododendron undergrowth on PRP 1, 2, and 5 (without intervention) was very vigorous with 81-97% cover percentage, a large number of monopodial woody stems amounting to 72,178-86,884 ha −1 , basal stem diameter of 4-82 mm, and varying mean height of individuals (141-224 cm). The coverage of rhododendrons on PRP 4 was 3% because its clearing was done 6 years ago, and no rhododendrons occur on PRP 3 due to clearing carried out 3 years ago. In terms of basal stem diameter structure, the first four diameter classes below 2.5 cm accounted for more than half (59.4%) of the total share of the stem number (Figure 4). In the left-shaped distribution of stem diameters, the largest representation of stems was in diameter class of 1.5-2 cm (13,800 stem ha −1 ). The share of rhododendrons with basal stem diameter >65 mm was only 0.4%.

Natural Regeneration
The number of recruits per hectare (natural regeneration) varied from 113 (PRP 1) to 105,075 (PRP 4); on the average eastern beech accounted for 97.0%, European hornbeam for 2.9%, and silver lime (Tilia tomentosa Moench.) for 0.1% ( Table 6). The cover percentage of rhododendron had a significant effect on the abundance of natural regeneration (F 1, 18 = 57.8, p < 0.001). The highest number of recruits occurred on the research plots with the zero or minimum cover percentage of rhododendron where this shrub was removed by clearing. Three years after clearing the recruit density on subplots of PRP 3 was 29,900-78,600 recruits ha −1 , and six years after clearing there were as many as 55,800-159,500 recruits ha −1 on subplots of PRP 4. A cover percentage of rhododendron abundance of natural regeneration on subplots of PRP 1, 2, and 5 was minimal i.e., 81-97% (100-1300 recruits ha −1 ).  Figure 5 shows a histogram of the height frequencies of natural regeneration of tree species on PRP 1-5. In case of PRP 3 (three years after rhododendron clearance), left-sided height distribution was observed with an exponential decline from a significant proportion of the first (h ≤ 20 cm) height class (12,088 recruits ha −1 ). After 6 years of rhododendron clearance (PRP 4), the most represented height class shifted to the right (40-60 cm; 10,100 recruits ha −1 ) because of dynamics. On other PRP (without rhododendron clearance), no height class reached the limit of 100 recruits ha −1 . On PRP 1 and 5, the height class distribution was uneven, in contrast to PRP 2 it was relatively balanced. The mean height of recruits was significantly taller (F (4, 7468) = 524.5, p < 0.001) on PRP 5 (221.5 cm ± 8.5 SE) and PRP 1 (196.7 cm ± 9.0 SE) than on PRP 2 (87.2 ± 5.0 SE) and PRP 4 (43.3 ± 0.4 SE); the smallest height was found on PRP 3 (23.2 ± 0.4 SE).
From a biodiversity perspective (Table 7), the spatial pattern of natural regeneration was significantly aggregated on all research plots (α = 1.986-21.053, R = 0.235-0.735, and ICS = 1.624-49.975) with the highest tendency of aggregation on PRP 3. PRP 5 is an exception where horizontal structure of recruits was insignificantly aggregated and/or random (according to α and ICS indices) due to their low number in distribution. It is also reported in the results obtained from the L-function ( Figure 6). Species richness, with monospecific or dominancy of eastern beech, evaluated by the D index of the relative measure of species diversity was zero on PRP 1, 2, and 5 and very low on PRP 3 and 4. It was the same in entropy H and species evenness indices. The vertical diversity was low on PRP 4 (A = 0.273) and medium on the other PRP (A = 0.362-0.487). The lowest values of TM h index, which describes the height differentiation, were obtained for PRP 4 (0.287) while the highest values were found on PRP 3 (0.651).   Table 5; * statistically significant (α = 0.05) for horizontal structure (α, R, and ICS index). Figure 7 illustrates relationships of the spatial pattern of recruits and the overstory on PRP 2 and 4 by means of the pair correlation function. These two plots were chosen where the highest number of recruits >1 m in height occurred and that represented both treatments: PRP 2-81% cover percentage of rhododendron (without intervention) and PRP 4-clearing the rhododendron. This figure indicates that the tree layer had a significantly positive influence (α = 0.05) on the spatial pattern of natural regeneration to a distance of 1.5 m on plots with the high cover percentage of rhododendron (PRP 2) where its vegetative reproduction (layering) was dominant. On the contrary, on plots where rhododendrons were removed, the tree layer had a negative influence on natural regeneration within 4 m (PRP 4).  Relationships of the spatial pattern of recruits and the overstory expressed by the pair correlation function on permanent research plots (PRP) 2 (rhododendron cleared) and permanent research plots 4 (no intervention); bold black line represents the pair cross function for real distances of trees; blue line on the level of g (r) = 1 represents the mean course for random spatial distribution of trees; and two thick blue curves represent 95% interval of reliability; if the observed value (black line) exceeds the upper (respectively lower) limit of the simulation interval, it indicates significant aggregation-positive relationship between two tested groups (respectively inhibition-negative relationship).

Relationships between Rhododendron, Tree Layer, and Natural Regeneration
Results of the PCA are presented in the form of an ordination diagram in Figure 8 and correlation matrix in Table 8. The first ordination axis explained 37.9% variability, the first two axes together explained 65.5% variability, and the first four axes together explained 86.2% variability in the data.
The first x-axis represented stand stocking (SDI) and canopy (crown closure and crown projection area) with height differentiation. The second y-axis represented total diversity index and mean height of natural regeneration with the number of recruits. Stand volume was significantly positively correlated with basal area (r = 0.97; p < 0.001) and mean stand height with DBH (r = 0.82; p < 0.001). Crown closure, SDI, and crown projection area increased with age, while these parameters were negatively correlated with height and diameter differentiation. The number of trees was significantly positively correlated with Arten-profile index (r = 0.52; p < 0.05) and total diversity index (r = 0.61; p < 0.01), while tree density was significantly negatively correlated with diameter differentiation (r = −0.58; p < 0.01). Total diversity also significantly decreased with increasing mean stand height (r = 0.63; p < 0.01) and DBH (r = 0.70; p < 0.001). Mean height of natural regeneration (recruits) was significantly positively correlated with the cover percentage of rhododendron (r = 0.67; p < 0.001), while these parameters were significantly negatively correlated with the density of natural regeneration (height: r = −0.52; p < 0.05; rhododendron: r = −0.84; p < 0.001). All these parameters were found independent (p > 0.05) of canopy (crown projection area, crown closure). Moreover, density of natural regeneration was significantly positively correlated with DBH (r = 0.57; p < 0.05) and significantly negatively with number of trees (r = −0.52; p < 0.05). Contribution of the mean stand age to explanation off variability in the data was relatively small. The difference within PRP was remarkable, especially for PRP 2, as marks of each record were relatively distant from each other, whereas marks for PRP 1 were fairly close to each other in the diagram. Research plots were significantly different from each other, as PRP 3 and 4 (with rhododendron cleared) with high density of natural regeneration occupied the upper part of the diagram, while higher stand density, rhododendron cover, vertical diversity, and mean height of natural regeneration were typical for PRP 1, 2, and 5 (lower part).   Table 3, Table 4, and Figure 8; * statistically significant at the level p < 0.05, ** p < 0.001.

Discussion
Eastern beech forest is one of the most important forest communities in the Black Sea Region [82] and Turkey, where it occupies 1.7 million hectares [14]. However, dynamics and structure of these valuable forest ecosystems are, in the long-term, significantly negatively affected by aggressive invasive pontic rhododendron [83,84]. A low species diversity was characteristic for the studied eastern beech stands in the BSR of Turkey in terms of species richness, heterogeneity, and evenness. Poor diversity is caused by strong competition of dense rhododendron cover and inappropriate management in its removal [85,86]. The horizontal structure of the tree layer was indistinctly aggregated or moderately regular while the random spatial pattern of trees across the research plots prevailed. Zenner et al. [87] also documented the prevailing random distribution of trees in near-natural stands of eastern beech in the Caspian Iran. The initial spatial pattern of near-natural beech stands is usually aggregated or randomly irregular and it turns to a moderately regular pattern in favorable conditions [88,89]. This may be due to gradual growth, competition among trees, and expansion of crowns as beech crowns are able to promptly react to the changes in the normal environmental condition [90].
The vertical structure of investigated stands varied from medium to very high. Diameter structure was also medium-diversified. Frequency distribution of diameter of the tree layer was mostly bell-shaped, which documents the stage of optimum [54]. Similarly, other studies show Gaussian distribution of diameters of eastern beech in the mature stands. Atalay [3] and Eşen et al. [4] also reported the relatively highest forest biodiversity in the oriental beech stands in the BSR within Turkish forests. A high biodiversity of near-natural forest stands with oriental beech was also reported in Iran [87,91]. These studies show that forests create the structural mosaic with alternation of particular developmental stages and phases. Sagheb-Talebi and Schütz [92] found irregular and diversified structure with a high proportion of large trees in the natural beech forests of the Caspian Iran.
The number of trees on the studied PRPs varies from 175 to 381 trees ha −1 . These numbers are roughly consistent with the data from the Elborz Mountain range along the southern Caspian Sea, i.e., 219-443 trees ha −1 [87]. In comparable stands of eastern beech in Iran, Akhavan et al. [93] reported 302-454 trees ha −1 and Sefidi et al. [17] reported 198 trees ha −1 on average. The highest stand volume was found on PRP 1, reaching 589 m 3 ha −1 , and the lowest stand volume was on PRP 2, i.e., 331 m 3 ha −1 . In eastern beech stands of identical growth conditions, Hassani and Amani [94] reported higher volumes (i.e., 562-678 m 3 ha −1 ) and Sefidi et al. [17] documented 386 m 3 ha −1 . Basal area of the studied stands was in the range of 28.8-46.1 m 2 ha −1 , and similarly Zenner et al. [87] found stand basal area varying from 33.7 to 44.1 m 2 ha −1 .
Natural regeneration of beech forests is strongly limited by dense undergrowth of aggressive invasive pontic rhododendron, which is very abundant in eastern beech stands in the entire western BSR [7,[95][96][97]. The number of monopodial woody stems varied from 72,178 to 86,884 stems ha −1 on our PRP, and this number was between the values reported by Eşen et al. [7], which were 48,000 and 91,000 stems ha −1 for R. ponticum and R. flavum, respectively. High rhododendron production was also observed from neighboring districts, where biomass storage capacity of pontic rhododendron reached up to 40,035 kg ha −1 with preference of sunny exposure and altitude to 500 m a. s. l. [98]. From this reason, pontic rhododendron can be (except positive effect of clearing on tree regeneration) used for fuel and in wood industry, especially for manufacture of density fiberboard panels [99,100]. Our study also reported dense rhododendron cover (81%-97%) and low tree regeneration (113-619 recruits ha −1 ). Only two research plots (PRP 3 and 4), where invasive rhododendron was cleared 3 and 6 years ago, respectively, contained a high number of tree regeneration (63,987-105,075 recruits ha −1 ). Comparable amount of beech natural regeneration (24,500 to 123,400 seedlings ha −1 ) was observed in the Sökü district in Western part of BSR [83]. Three to four individuals per m 2 may be considered as a sufficient number of tree regeneration for the prosperous development of stands [101]. From this perspective, only tree regeneration on PRP 3 and 4 with 0-3% cover percentage of rhododendron may be considered as sufficient. On the average, eastern beech accounted for 97% of natural regeneration (recruits) on the studied research plots. Because site preparation was done in seed mast year of beech, the seed status of other tree species was not taken into consideration. As a result, regeneration was purely dominated by beech seedlings, but other tree species consisted of admixture of mature stands. However, rhododendron had an opposite influence on the mean height of recruits. Due to competitiveness of this woody shrub, the height of recruits was significantly taller on the research plots with a high cover percentage of rhododendron (87.2-221.5 cm) compared to the recruit heights (23.2-43.3 cm) on PRPs where rhododendron was removed. While the influence of rhododendron cover (light competition) on regeneration was strong, influence of overstory canopy closure on recruit number and height was not proved.
The horizontal structure of natural regeneration was distinctly aggregated on all studied PRPs. A similar pattern was also found in comparable oriental beech stands [93]. The situation is identical in European beech (Fagus sylvatica L.) [19,102]. A tendency towards aggregation of seedlings and in the lowerstory is a result of prosperous growth of natural regeneration in canopy gaps that originated after dead trees [23,103]. As for the spatial pattern (contrary to recruit number and height), parent trees had a significant negative influence on the individual recruits within a 4 m radius from tree stem on the PRPs, where rhododendron cover with generative regeneration prevailed. Vacek et al. [19] reported that distribution of parent trees negatively affected the placement of natural regeneration of European beech within a 2.5 m radius. This negative spatial interaction between natural regeneration and tree layer in small spacing may be caused by insufficient light conditions for growth of seedlings closer to mature tree. However, trees had positive influence on natural regeneration to the distance of 1.5 m on the research plots with dense rhododendron cover. This is caused by rich root sprouting of eastern beech [17]. Natural regeneration, especially from the seed without support from the parent plant, has no chance to survive under the dense rhododendron cover compared to the sprouts or suckers [7]. The species diversity of natural regeneration was low on the study site; however, it was higher on PRP 3 and 4, where rhododendron was removed. On other plots, no other tree species (except beech) occurred. Removal of rhododendron may give an opportunity to the other tree species to establish regeneration. However, during the rhododendron clearance, the other tree species are usually removed from the beech forests, and frequent use of bulldozers had a negative effect on the seed bank during the soil removal [86].
The mineral soil exposure during site preparation is often recommended for sufficient natural regeneration because the litter layer does not provide a favorable seedbed for beech seeds [96,104]. In this context, Agestam et al. [105] stated that the suitable preparation of site conditions, and especially of the soil surface with the exposure of mineral contents, increased the germination rate of beech seeds compared to sites without soil preparation in southern Sweden. One year after the bulldozer removal of rhododendron undergrowth along with thickness of 5-10 cm soil layer in beech forests of the BSR, Yildiz et al. [95] and Sarginci [106] reported 300,000-700,000 beech seedlings ha −1 . Such a method of soil preparation is more suitable from short-time economic aspects than manual cutting or clearance of rhododendron. At the same time, it leads to the accelerated mineralization causing nutrient loss, decrease in microbial activity, and soil compaction [95,101,107]. Soil compaction can also hinder the root growth, decrease tolerance to drought, and limit the uptake of immobile soil nutrients [108]. For these reasons, and especially from the perspective of the sustainable forest management, manual cutting and clearance of rhododendron undergrowth seems better than its removal by a bulldozer with great disturbance of the soil surface. In addition, during site preparation, scraping the forest floor by a rake-equipped bulldozer significantly disturbed the seed banks. A similar opinion was published by [101], who supported manual clearance of rhododendron and its burning because of their efficiency for the long-term productivity of oriental beech stands in Turkey.

Conclusions
The findings of this study document an invasive spread of pontic rhododendron as a limiting factor for prosperous regeneration in beech forests in the Black Sea Region, which can be further enhanced in future due to climate change. This aggressive invasive shrub has a long-term significant negative effect on natural regeneration and subsequent dynamics and structure of forest ecosystems. Rhododendron cover significantly affected spatial patterns of natural regeneration, considering its mean height and especially density. Number of tree regeneration was 309× (188×) higher in forest stands after 6 (3) years after rhododendron clearance compared to stands without intervention with a mean annual increase in the regeneration number by 19.4 thousand recruits ha −1 y −1 . On the other hand, tree canopy had no significant effect on density of tree regeneration and rhododendron cover. The natural regeneration of these forests will be very problematic without the efficient clearance of rhododendron by effective methods minimally disturbing the soil environment and seed bank, while blockage of the successional stages of rhododendron may occur during disintegration of the beech stands. The appropriate methods of rhododendron removal and the regeneration techniques have significant positive effect on continuity of forest dynamics and thus protect the threatened beech forest ecosystems, but it can also increase forest species diversity. Other tree species (except eastern beech) occurred in natural regeneration only in stands after rhododendron clearance. Therefore, manual clearance of rhododendron (with the possibility of use in wood industry) at the beginning of the regeneration period is recommended in managed forests.