Disentangling the Ecological Determinants of Species and Functional Trait Diversity in Herb-Layer Plant Communities in European Temperate Forests

: Forest herb-layer vegetation responds sensitively to environmental conditions. This paper compares drivers of both taxonomic, i.e., species richness, cover and evenness, and functional herb-layer diversity, i.e., the diversity of clonal, bud bank and leaf-height-seed plant traits. We investigated the dependence of herb-layer diversity on ecological determinants related to soil properties, climatic parameters, forest stand characteristics, and topographic and abiotic and biotic factors associated with forest ﬂoor structure. The study was conducted in different forest types in Slovenia, using vegetation and environmental data from 50 monitoring plots (400 m 2 each) belonging to the ICP Forests Level I and II network. The main objective was to ﬁrst identify signiﬁcant ecological predictors and then quantify their relative importance. Species richness was strongly determined by forest stand characteristics, such as richness of the shrub layer, tree layer shade-casting ability as a proxy for light availability and tree species composition. It showed a clear positive relation to soil pH. Variation in herb-layer cover was also best explained by forest stand characteristics and, to a lesser extent, by structural factors such as moss cover. Species evenness was associated with tree species composition, shrub layer cover and soil pH. Various ecological determinants were decisive for the diversity of below-ground traits, i.e., clonal and bud bank traits. For these two trait groups we observed a substantial climatic signal that was completely absent for taxonomy-based measures of diversity. In contrast, above-ground leaf-height-seed (LHS) traits were driven exclusively by soil reaction and nitrogen availability. In synthesis, local stand characteristics and soil properties acted as the main controlling factors for both species and trait diversity in herb-layer communities across Slovenia, conﬁrming many previous studies. Our ﬁndings suggest that the taxonomic and functional facets of herb-layer vegetation are mainly inﬂuenced by a similar set of ecological determinants. However, their relative importance varies among individual taxonomy- and functional trait-based diversity measures. Integrating multi-faceted approaches can provide complementary information on patterns of herb-layer diversity in European forest plant communities.


Introduction
Examining plant-environment linkages is a longstanding quest in vegetation ecology [1,2]. Despite the relatively well-established consensus among forest ecologists that local forest site factors should play a critical role in driving herb-layer vegetation (composed of both herbaceous plant species and woody seedlings and saplings) diversity, uncertainties remain whether the same external or internal features consistently affect different facets of species organization [3]. Plant community diversity and composition are influenced by ecological determinants operating over different spatial scales [4]. At broader scales, climate is often considered one of the main explanatory variables for the diversity of forest plant communities [5,6]. At finer spatial scales, however, other interacting abiotic conditions and biotic attributes, such as forest stand characteristics (e.g., tree layer composition and diversity; [7,8]) and soil properties (e.g., physical and chemical properties of soil; [9,10]) likely come into play in shaping patterns of floristic and trait diversity. For example, the authors of [6] demonstrated that forest understory vegetation is more related to soil than to climate towards the cold distribution margin of European beech. Furthermore, topography (aspect and slope; [11,12]) and associated microclimatic conditions (air temperature and humidity; [13,14]) can also be decisive for the diversity and abundance of forest vegetation. In comparison to non-forest vegetation (e.g., grasslands) with less vertically stratified structure, the interplay between different vegetation layers in the forest ecosystem (i.e., tree, shrub and herb layer) additionally increases the complexity of the dependencies of herb-layer vegetation on ecological factors coming from above (light) and below (soil moisture and nutrients). These determinants are particularly important for the abundance and diversity of herbaceous species [15][16][17], whereas the woody component of the herb layer depends to a greater extent on factors such as overstory seed sources or management preferences [18,19].
The forest understory can be seen as a resource-limited habitat [20]. Herbaceous and woody plants underneath overstory canopy are influenced by overstory composition and structure through modifications of essential resource (light, water and soil nutrients) availability and heterogeneity [7,21], but also by other effects, such as the characteristics of the litter layer [4,22]. The filtering of resources induced by the overstory tree layer directly impacts herb-layer diversity. One classic example is the prevention of the establishment of more light-demanding species due to the low-light understory environment in closed canopy mature stands composed of late-successional tree species [15,23]. Canopy characteristics regulate light availability at the forest floor [24]. The effects of the aforementioned factors on herb-layer assemblages are often interrelated [21], further complicating our capacity to assess the contribution of each environmental variable to taxonomic or functional diversity. For instance, forest soil properties depend on the dominant tree species [25] and macroclimate [6].
Comparisons of different forest types with differential abiotic and biotic conditions offer a good opportunity to study ecological factors controlling the diversity of forest herb-layer vegetation. However, within the pan-European programme for monitoring forest health (i.e., the International Co-operative Programme on the Assessment and Monitoring of Air Pollution Effects on Forests-ICP Forests; [26][27][28]), previous studies rarely focused on ground vegetation. The herb layer should be a central part of forest diversity studies because herbaceous species majorly contribute to vascular plant species richness and overall forest biodiversity [29][30][31] and play a pervasive role in ecosystem functioning [21,32]. Vegetation records from ICP Forests monitoring plots (Level I or Level II; [33]) have been more frequently used for the research of possibilities and accuracy of different survey methods [34,35] or impacts of atmospheric nitrogen deposition [36][37][38] and forest disturbances [39]. Only recently have authors [40,41] addressed variation in forest understories along environmental gradients from a trait-based perspective.
In the last decade, studies dealing with the ecological determinants of forest herb-layer diversity have expressed the need for the implementation of multi-faceted approaches, where the taxonomic aspect (species identity) should be complemented at least with a functional diversity facet [42,43] or the phylogenetic structure of vegetation [44,45]. This is relevant for two reasons: (1) plant functional traits offer valuable insights into the way a plant is linked through its morphology and physiology to the environment and its resources, and how it influences ecosystem functioning and the provision of ecosystem services [19,32,46,47], and (2) abiotic and biotic drivers of species-level diversity may differ from those that significantly affect trait diversity. However, simultaneous investigation of species diversity and interspecific trait dissimilarity along environmental gradients [47,48] has been scarcely utilized to the extent that would give informative results or even have applied value. In addition, trait-based studies are largely over-represented by traits con-nected to plant performance through above-ground processes (e.g., competition for light capture, dispersal ability). This limits our understanding of general species assemblage organization and explicitly calls for stronger integration of below-ground traits into the functional traits approach in plant ecology [49], as they add a new, independent dimension to the plant economic spectrum [50,51]. Traits related to different functions should show different responses along particular ecological gradients, but our understanding still lacks generalizable conclusions.
To fill this knowledge gap, we analyzed a comprehensive dataset from ICP Forests monitoring plots across Slovenian forest ecosystems. Our study has two main goals: (1) to assess the relationship between different measures related to taxonomy-based diversity (herb-layer species richness, cover and evenness) and groups of functional traits (clonal traits, bud bank traits and leaf-height-seed trait scheme) along various forest soil, climate, forest stand and topographic and structural factors and (2) to test the relative importance of these explanatory variables for different taxonomy-and trait-based measures of diversity. We hypothesized that a different set of ecological variables would be linked to taxonomic diversity compared to functional diversity, meaning that these diversity facets are likely best explained by a different suite of ecological factors.

Study Area, Vegetation Sampling and Data Collection
The study area encompasses 39 monitoring plots of ICP Forests Level I and 11 monitoring plots belonging to ICP Forests Level II (Intensive Monitoring Programme), i.e., 50 plots in total ( Figure 1). ICP Forests is a transnational network for monitoring forest health status in Europe [27]. Level I plots are systematically distributed in a 16 km × 16 km grid across Slovenia [52]. Level II plots were established in 11 different locations, representative of the heterogeneity of Slovenian forest types [53]. Detailed information on Intensive Monitoring plots is reported in the work of [39] and in [54]. Collectively, Level I and Level II plots were located in forest sites with a wide variety of climatic, geological, edaphic and topographical conditions as well as a broad range of forest vegetation composition. The altitude of our plots ranged from 160 to 1490 m, mean annual temperature from 3.2 to 11.7 • C and mean annual precipitation from 791 to 2499 mm (the climatic parameters were taken from the WorldClim database, https://www.worldclim.org/ [55], accessed on 15 November 2020). Forest ecosystems cover 58% of Slovenia's land area, which makes it the third most forested country in Europe [56,57]. Due to its specific geographic position in the transition between the major geographical units of South-eastern and Central Europe (the Alps, the Dinaric Mountains, the Mediterranean Region, the Pannonian Basin) and diverse climatic and topographic conditions, Slovenia exhibits a wealth of forest plant communities [58]. A diverse mosaic of forest types includes lowland floodplain forests, widespread mixed forests at intermediate altitudes and mountainous coniferous forests at higher altitudes [39,59]. The majority (76%) of forests in Slovenia are private property, 24% of forests are public. Fundamental principles of forest treatment and management include sustainability, close-to-nature management and multi-purpose management.
Forest vegetation in monitoring plots was sampled in the period 2006-2007 (Level I plots) and in 2009 (Level II plots), following harmonized ICP Forests protocols [33]. Surveys were made during summer months (June, July) when all forest vegetation layers were fully developed. The size of the sampling area was 400 m 2 in all plots. All vascular plant species in the tree (woody individuals taller than 5 m), shrub (woody individuals from 0.5 to 5 m tall) and herb layer (all herbaceous species and woody individuals lower than 0.5 m) were recorded. The tree layer contains trees and shrubs, and the herb layer contains herbaceous (forbs, graminoids, ferns) and woody species. The abundance of encountered species in each layer was visually estimated using a modified Barkman scale [60]. Species nomenclature followed the work of [61][62][63]. During the vegetation survey, the following parameters were also estimated: the overall cover (%) of each vegetation layer (tree, shrub, herb and moss layer-mosses and other bryophytes covering the forest soil), cover (%) of outcropping rock and rock fragments and cover (%) of woody debris. For each plot, topographic information was recorded, such as altitude (m above sea level), aspect (azimuth degrees) and slope (vertical degrees).

Soil Properties
The soil data used for the present investigation were collected during the BioSoil biodiversity project [64]. We made pedological profiles in the soil, sampled them by probing and described their morphological characteristics. In each monitoring plot, soils were sampled at five sampling points, which were then pooled to obtain aggregated samples [65]. At each sampling point, soil samples were collected from the organic layer (Oh, Ol, Of) and different depths of mineral soils: 0-5 cm layer, 5-10 cm layer, 10-20 cm layer, 20-40 cm layer and 20-80 cm layer. The following variables were measured for soil samples from different layers: pH (CaCl 2 ), total nitrogen content (g/kg), soil texture, soil organic carbon (g/kg), moisture content (%), cation exchange capacity (cmoL/kg), exchangeable calcium (cmoL/kg), bulk density (kg/m 3 ), extractable phosphorus (mg/kg) and many others. Soil variables were measured according to standard procedures [66] and ICP Forests protocols [53,67].
Since our complete dataset contained a few tens of different soil variables, we investigated (inter)correlations between these variables to detect variables containing similar information. Based on this, we excluded highly correlated variables and selected five different soil variables related to soil pH (pH of the organic layer, pH of the upper mineral soil-depth up to 20 cm), soil productivity (total nitrogen in the organic layer, total nitrogen in the upper mineral soil-depth up to 20 cm) and soil texture (proportion (%) of clay in the first 20 cm of mineral soil). For pH, nitrogen availability and clay content, the values from different layers were aggregated and averaged: Oh, Ol and Of for the organic layer and three layers (depth 0-5 cm, depth 5-10 cm and depth 10-20 cm) for the upper mineral soil. The organic and upper mineral soil layers were included because topsoil is considered to be most important for plant species in the herb layer [68].

Climatic Parameters
Climatic conditions for each plot were compiled from the WorldClim global database, version 2.1 [69], with a spatial resolution of 30 s (~1 km 2 ). Among 19 standard bioclimatic variables, two variables related to temperature (i.e., temperature seasonality-BIO4, maximum temperature of the warmest month-BIO5) and two variables associated with precipitation (i.e., precipitation seasonality-BIO15, precipitation of the wettest month-BIO13) were chosen. We additionally included data on total potential evapotranspiration (PET, representing atmospheric water demand), which was sourced from the Global Aridity and PET Database [70]. Overall, five different climatic parameters were selected. The raster package in R [71] was used for downloading the data from both open-source databases.

Forest Stand Characteristics
We collected data on several stand characteristics related to the tree and shrub layer. From vegetation surveys [33,66], tree layer richness (number of different tree species in the tree layer) and the proportion of broadleaf species (proportion of broadleaves in the total cover of all species in the tree layer) were calculated. Furthermore, shade-casting ability (SCA; [23,72]) was calculated as a community-weighted mean [73] where shade production of tree species [19] was used and weighted by the relative abundance of species in the tree layer. The authors of [23] suggested that accounting for compositional canopy characteristics could improve predictions for the light-demand signature of the forest understories. Vegetation data were also used to calculate the richness of the shrub layer (number of different woody species in the shrub layer) and estimated values of shrub layer cover. Overall, five different plot-level stand characteristics were selected (Table 1).

Topographic and Structural Factors
In this group of explanatory variables, different abiotic and biotic variables associated with local topography and forest floor structural features were considered, namely altitude, slope aspect, slope angle, cover of rocks on the surface, cover of woody debris and cover of the moss layer. Topographic factors (aspect and slope) were used to calculate the heat load index (H i ) with the following equation [74]: where aspect is in degrees azimuth and slope is in degrees. Overall, five different plot-level topographic and structural factors were included in the analysis. A detailed description of 20 numerical predictors, divided into four groups, is given in Table 1.

Response Variables
All response variables were derived from the herb-layer data. We distinguished between three taxonomy-based measures of diversity: herb-layer species richness, herblayer cover (%) and species evenness. Species richness means the number of different plant species and is a key variable in measuring the diversity of ecological communities [35]. Herb-layer cover was visually estimated during vegetation sampling. Species evenness was calculated using Pielou's equation [75]: The diversity function in the R package vegan [76] was implemented for this purpose. The Shannon or Shannon-Weaver diversity index is defined as H = −Σ i p i ln p i , where p i is the proportional abundance of species i and ln is the natural logarithm. Species evenness is a description of the distribution of relative abundance across the species in a community, i.e., abundance equality of species [77]. For the functional trait-based measures of diversity, three different groups of plant traits were considered ( Table 2). We selected a suite of functional traits involved in both above-and below-ground processes. First, we gathered data on clonal traits: clonality, clonal index, and the type and persistence of the clonal growth organ. Second, two bud bank traits were included: bud bank size and bud bank depth. Traits associated with clonality and bud banks provide information on space occupancy and vegetative resprouting [41], both affecting plant local persistence as shown in the work of [78]. Clonal and bud bank traits also pose a prominent signature in the context of species responses to different types of disturbance [79]. Lastly, traits connected to plant resource acquisition and the rate of their investment (specific leaf area-SLA), competitive ability (plant height-H) and dispersal capacity (seed mass-SM) represented the third group of functional traits. These are commonly described by the leaf-height-seed (LHS) strategy scheme [80,81]. The authors of [82] illustrated that intraspecific variability of SLA fosters the persistence of forest specialists across a light availability gradient in beech understories. Clonal and bud bank traits were extracted from the CLO-PLA database [83,84]. Useful definitions for these traits are available in the Pladias online database [85] (https://pladias. cz/en/ accessed on 15 December 2020). Values of SLA, H and SM were obtained from the LEDA database [86]. For H, only herbaceous species were included in the analyses as plant height for woody (tree) species is parameterized for adult trees. Seed mass for pteridophytes (sexual reproduction via spores) was set to zero, similarly to [4].
We used Rao's quadratic entropy (RaoQ; [87]) as a measure of functional diversity, i.e., how different coexisting plant species are in terms of their traits [88]. RaoQ is an index of functional dissimilarity, expressed as the sum of the abundance-weighted pairwise differences between species in a community [73,87]. It was calculated separately for the three trait groups, i.e., trait combinations: clonal traits, bud bank traits and LHS traits. As the group of clonal traits consisted of both qualitative and quantitative variables, the procedure described by [87] was considered, i.e., using Gower distance [89] with the analytical modification proposed by [90]. The FD package [91] was employed for RaoQ computation (function dbFD). Some numerical traits were log-transformed prior to calculation as their values ranged between several orders of magnitude and to improve the distribution of residuals. Plot-level values of the RaoQ of each trait combination were used as a response variable in further data analyses.
We checked for interdependence between response variables by calculating the Pearson correlation coefficient (Supplementary Materials, Table S1).

Statistical Analyses
All statistical analyses and visualizations were performed in the R programming language for statistical computing and graphics, version 3.5.2 [92].
Generalized additive models (GAMs; [93]), a non-parametric extension of generalized linear models [94], were used to identify the most significant predictors for each herb-layer response variable: species richness, cover, evenness, RaoQ for clonal traits, RaoQ for bud bank traits and RaoQ for LHS traits. Before the stat checked.
Given that our complete dataset contained 50 plots (observations) and 20 predictors (explanatory variables), we avoided using all predictors at the same time. We thus prereduced the number of predictors to be included in the GAM. Modelling was performed in two steps, using the gam function in the mgcv R package [95]. We first fitted "local" GAMs with predictors from the same group (i.e., soil parameters, climatic parameters, stand characteristics, topographic and structural factors). This resulted in four different local models for each response variable, 24 local models in total (Supplementary Materials, Figure S1). The number of potential predictors in each group was balanced (i.e., five per group).
Below is an example of the R code for local GAM fitting: where y represents the response variable, while arguments from var 1 to var 5 are covariates, i.e., predictor variables. To avoid multicollinearity, GAM model selection was undertaken by restricted maximum likelihood (REML) with automatic variable selection based on adding penalties to successive variables [96]. For a detailed explanation of the specified arguments in the code above, see the mgcv package manual file available on the CRAN Repository (https://cran.r-project.org/, accessed on 15 January 2021). Based on these local GAMs, predictors from different groups were ranked in terms of their overall effect size (F-statistics) and statistical power, from most to least significant. In the second modelling step, we constructed a "global" (final) GAM for each response variable, where we selected the five most significant predictors based upon their ranking derived from the local models. In terms of model specifications, global GAMs (one for each response variable) were fitted in the exact same way as local models (formula above). In all cases, the function gam.check in the mgcv package was used for model diagnostics with appropriate adjustment of the number of nods in the smoothed functions when necessary. The adjusted coefficient of determination (adj. R 2 ) was used as an overall measure of goodness-of-fit for the final GAMs. The p-value in GAM analysis for each explanatory variable indicates whether its effect within the model is significantly different than zero.
Before model construction, explanatory variables were inspected for their distribution. When significant departures from normality (checked with histograms and the Shapiro-Wilk test) were detected, appropriate data transformations were performed. For this we used the R package bestNormalize [97]. To further improve the interpretability of explanatory variables in terms of their relative importance, standardization of predictors was implemented, resulting in scaled and centered values with mean = 0 and standard deviation = 1. Response variables were also checked for their distribution. No significant deviation from normality was detected. Therefore, Gaussian error distribution and the identical link function in model specification (see formula above) was an appropriate choice.
Variance partitioning [98] was used to identify the relative role of groups of ecological factors on the response variables. This statistical approach has frequently been used to investigate how different groups of ecological factors affect species and trait distributions [99]. Only predictors included in the final GAMs were used. The function varpart in the vegan R package [76] served for variance partitioning analysis.

Herb-Layer Vegetation Status across Monitoring Plots
We recorded a total of 402 different vascular plant species in the herb layer, of which 80 were woody species (trees, shrubs, woody climbers) and 322 were herbaceous plants (forbs, graminoids, ferns). The mean (± standard deviation) plot-level species richness of the herb layer was 44.0 ± 20.5. The five most frequent plant species recorded in the herb layer were Acer pseudoplatanus (present in 76% of all plots), Picea abies (64%), Fagus sylvatica (62%), Mycelis muralis (60%) and Oxalis acetosella (56%). The herb-layer communities across 50 monitoring plots exhibited a profound range of both taxonomy-and trait-based diversity measures (Table 3). Table 3. Descriptive statistics for three taxonomy-based measures of diversity (species richness, cover, evenness) and three variables related to functional diversity (clonal, bud bank and leaf-height-seed (LHS) traits), based on 50 herb-layer forest plant communities. For traits, Rao's quadratic entropy (RaoQ) was used as a measure of functional diversity.

Herb-Layer Species Richness
For herb-layer species richness, all five predictors included in the final GAM were significant (Figure 2) with the highest explanatory power (adjusted R 2 = 65.1%) among all final models. Three variables related to stand characteristics (shrub layer richness, tree-layer SCA and the proportion of broadleaves in the tree layer cover) and one related to soil properties (pH of organic soil) and structural factors (cover of rocks on the surface) were selected. Herb-layer species richness increased along the gradient of shrub layer richness, organic soil pH and surface rockiness. A similar pattern was observed for the gradient of tree-layer SCA but with an evident plateau at higher values of tree-layer SCA. In contrast, herb-layer species richness decreased with the proportion of broadleaves in the tree layer cover (Figure 2).

Herb-Layer Cover
For the total herb-layer cover, three predictors (shrub-layer species richness, moss cover, proportion of broadleaves in the tree layer) were statistically significant, whereas the other two (temperature seasonality-BIO4 and tree-layer SCA) were marginally significant ( Figure 3). Herb-layer cover increased with shrub layer richness. The relationship between herb-layer cover and moss cover was hump-shaped with the lowest values of herb-layer cover in the middle of the gradient. A similar pattern was observed for temperature seasonality and tree-layer SCA. Herb-layer cover decreased along the gradient of the proportion of broadleaves in the tree layer ( Figure 3). The final GAM explained 32.7% of the variation in herb-layer cover.

Herb-Layer Species Evenness
For herb-layer evenness, the most significant predictor was the proportion of broadleaves in the tree layer, followed by the total cover of the shrub layer and organic soil pH (Figure 4). Total nitrogen content in organic soil was marginally significant. The highest herb-layer evenness occurred in the middle of the gradient for the proportion of broadleaves in the tree layer (unimodal curve). It increased with total shrub layer cover and with organic soil pH but decreased along the gradient of nitrogen availability in organic soil (Figure 4). The final GAM explained 29.9% of the variation in herb-layer evenness.  Table 1.

Diversity of Clonal Traits
For herb-layer diversity of clonal traits, the most significant predictor was shrublayer species richness, followed by precipitation seasonality (BIO15) and organic soil pH ( Figure 5). The content of clay in the mineral soil was marginally significant. The highest diversity of clonal traits was observed for plots with an intermediate level of shrub layer richness. It increased with precipitation seasonality and also with organic soil pH. However, trait diversity decreased along the gradient of the clay content in the mineral soil ( Figure 5). The final GAM explained 43.6% of the clonal trait variation. The estimated relationships between the herb-layer cover and five explanatory variables included in the final GAM: shrub layer richness, moss layer cover, proportion of broadleaves in the tree layer, temperature seasonality (BIO4) and tree-layer shade-casting ability (SCA). The order of predictors reflects their effect size and significance: * p < 0.05, m (marginally significant)-p < 0.1. The confidence interval (95%) is indicated with the light-blue ribbon around each curve. For additional definition of predictors, see Table 1.

Diversity of Bud Bank Traits
In case of herb-layer bud bank traits, one ecological predictor from each group was selected in the final model: altitude (topographic factor), pH of the organic soil layer (soil attribute), tree-layer SCA (stand characteristic) and precipitation seasonality-BIO15 (climatic parameter) ( Figure 6). The diversity of bud bank traits decreased with altitude and tree-layer SCA while exhibiting an increase with the pH of the organic soil layer. The relationship between bud bank trait diversity and precipitation seasonality was curvilinear, with the lowest functional diversity in the middle of the gradient (Figure 6). The final GAM explained 34.9% of the bud bank trait variation. . The estimated relationships between the herb-layer species evenness and four explanatory variables included in the final GAM: proportion of broadleaves in the tree layer, shrub layer cover, pH of the organic soil layer and total nitrogen content in the organic soil layer. The order of predictors reflects their effect size and significance: ** p < 0.01, * p < 0.05, m (marginally significant)-p < 0.1. The confidence interval (95%) is indicated with the light-blue ribbon around each curve. For additional definition of predictors, see Table 1.

Diversity of Leaf-Height-Seed Traits
For herb-layer LHS trait diversity, only two predictors were significant in the final model, both related to soil properties: pH of the organic soil layer and nitrogen content in the organic soil layer (Figure 7). The diversity of LHS traits showed a variable pattern along the organic soil pH gradient. It peaked at rather low values of pH and then generally decreased along the gradient. Herb-layer communities with a higher diversity of LHS traits were found on plots with lower nitrogen content in the organic soil layer, whereas on more productive soil, LHS trait diversity was generally lower (Figure 7). Despite a small number of significant predictors in the final GAM, the explanatory power of the model was relatively high (adjusted R 2 = 39.4%). Figure 5. The estimated relationships between the diversity of clonal traits and four explanatory variables included in the final GAM: shrub layer richness, precipitation seasonality (BIO15), pH of the organic soil layer and clay content in the upper mineral soil. The order of predictors reflects their effect size and significance: *** p < 0.001, * p < 0.05, m (marginally significant)-p < 0.1. The confidence interval (95%) is indicated with the light-blue ribbon around each curve. For additional definition of predictors, see Table 1.

Variance Partitioning
The analysis of variance partitioning (Figure 8) revealed that the relative contribution (expressed in %) of soil properties towards explaining the variation in response variables was as follows: clonal traits (10.4%) > LHS traits (9.8%) > species richness (6.3%) > evenness (5.7%) > bud bank traits (5.5%) > herb-layer cover (0.0%). Climatic parameters explained 6.6% of the variation in clonal traits, while all other GAMs did not include any of the climatic parameters. In contrast, stand characteristics explained a substantial proportion of variation in species richness (25.0%) > clonal traits (9.4%) > herb-layer cover (8.3%) > evenness (6.6%) > bud bank traits (6.0%), but did not contribute to explaining the variation in LHS traits (0.0%). Topographic and structural factors proved to have an important role in explaining the variation in bud bank traits (9.0%), species richness (3.2%) and herblayer cover (1.2%). However, these ecological drivers were not included in the GAMs for evenness, clonal traits and LHS traits (Figure 8). Figure 6. The estimated relationships between the diversity of bud bank traits and four explanatory variables included in the final GAM: altitude, pH of the organic soil layer, tree-layer shade-casting ability (SCA) and precipitation seasonality (BIO15). The order of predictors reflects their effect size and significance: ** p < 0.01, * p < 0.05. The confidence interval (95%) is indicated with the light-blue ribbon around each curve. For additional definition of predictors, see Table 1.
Averaging these results together, forest stand characteristics on average explained the highest fraction (9.3%) of variation in all six response variables considered. The second most explanatory group of predictors was soil properties (6.3%), followed by topographic and structural factors (2.2%) and climatic parameters (1.1%). Looking at the three taxonomy-based diversity measures separately, the order of predictor groups from most to least important was as follows: stand characteristics (13.5%) > soil properties (4.0%) > topographic and structural factors (1.5%). Climatic parameters did not contribute to explaining the taxonomy-based measures. In the case of the three response variables describing functional trait diversity, the order changed in favor of soil properties (8.6%). The variation in trait groups was also strongly affected by stand characteristics (5.1%) and topographic and structural factors (3.0%) but less so by climatic parameters (2.2%).

Species-Environment Relationships
Given that our study addresses the majority of forest types found in Central and SE Europe, the results are meaningful for similar European forests and also temperate forests in, e.g., Asia [100,101], North America [102][103][104] and other parts of the world [105]. We distinguished between three taxonomic response variables of herb-layer vegetation: species richness, herb-layer cover and species evenness. The first two have frequently been explored in many different forest types in the temperate zone around the globe [22,101,106], whereas species evenness has been examined to a much lesser degree.
Our findings indicate that shrub layer characteristics are important determinants for all three taxonomic response variables. Shrub layer species richness proved to be the most effective predictor of herb-layer richness and cover. This forest stand characteristic showed a positive effect on both variables. However, the increase in herb-layer richness along the gradient of shrub layer richness hit a plateau at high shrub-layer richness, suggesting competitive interactions between these two forest strata. Shrub-layer cover also proved significant for herb-layer species evenness. This again indicates the regulative role of the shrub layer, i.e., higher shrub cover prevents the domination of one or a few herblayer species. In temperate forests, however, the shrub layer is often overlooked as it is assumed that shrubs are rather sparsely distributed under a dense tree canopy layer. This might be true in most cases, but the shrub layer likely acts as a direct filter for the herb-layer species underneath, particularly influencing light availability and the forest microclimate [107]. Although the shrub layer reduces forest floor light availability and contributes to higher litter production, its presence or absence may create highly variable microenvironments [108,109].
In contrast to the relative importance of the shrub layer in our monitoring plots, previous studies suggested that forest herb-layer vegetation is largely controlled by treelayer characteristics. The authors of [106] found that tree species richness acts as a positive factor for herb-layer species richness and cover in mixed temperate forests in Hungary. None of our final GAMs included tree species richness as an explanatory variable. Instead, the proportion of broadleaves in the tree layer (a proxy for tree species composition) was a significant predictor in all three GAMs for taxonomy-based measures of diversity. The results indicated that a higher proportion of broadleaves tends to decrease herb-layer species richness and cover. At first glance, such patterns run counter to prevailing evidence in forest ecology.
In temperate and boreal zones, coniferous forests generally provide less diversified vascular understories than broadleaved forests [7]. This incongruence can be somewhat explained by the fact that our monitoring plots included some sites with pure conifer stands (e.g., Picea abies) growing on mixed parent material, which contributed to the coexistence of acidophilous species (present due to the thick litter layer from spruce trees that lowered the pH of the topsoil) and basophilous plants (occurring on less acidic microsites). On the other extreme, some stands were composed only of broadleaved tree species but had relatively low herb-layer richness. One typical example is the phytosociological association of Blechno-Fagetum, a species-poor beech community due to soils with low pH values and impoverished nutrients in the upper horizons. The observed importance of broadleaves (deciduous trees) might be related to seasonal variations of overstory canopy cover and associated herb-layer development. The compositional differences between spring and summer are relevant. It is important to note, however, that the summer aspect of herb-layer vegetation (capturing most of the seasonal variability of ground flora in temperate forests) was considered here. This sampling included both spring and summer species and is representative for these types of temperate forest.
Tree species composition was the most expressed predictor for herb-layer species evenness. The relationship between these two variables was evidently unimodal, suggesting that mixed forest stands support the highest herb-layer evenness. Mixed overstory canopy comprised of both broadleaves and conifers exhibits higher spatial and temporal variation of resources in the understory layer [21]. In contrast, the majority of the monodominant overstories in our dataset were characterized by the dominance of a single or a few herb-layer species.
Another important ecological determinant for herb-layer diversity was tree layer shade-casting ability, a proxy for light conditions at the forest floor. However, the results of our study suggest that species richness increased along the tree-layer shade-casting ability gradient, at least to a certain point-at the highest tree-layer SCA values, the non-linear curve started to bend down. This observation is slightly counterintuitive to the general expectation that more light at the forest floor would result in higher plant richness [110,111]. Light is the most important limiting factor in the majority of forest ecosystems [112], and shade tolerance affects the ability of plants to cope with other environmental stressors [20].
However, the effect of light is not straightforward because the diversity of the herb layer may not necessarily be in correlation with light availability. The authors of [13] noted that photosynthetically active radiation was of secondary importance for the distribution of widespread temperate forest herbs in beech forests on calcareous soils. In the study of [113], light conditions were the most important drivers in acidophilic beech-oak forests, while on more neutral forest sites different soil properties proved more determinant. Other studies (e.g., [25]) have also concluded that the effects of light availability depend on edaphic conditions.
One possible explanation is that on certain sites (pine forests in our case), higher light availability supports the expansion of plant species adapted to a high-light environment (e.g., grasses). Their dominance can lead to the competitive exclusion of less adapted and competitively inferior species, resulting in the impoverishment of community-level richness. Another reason for this pattern may, at least in part, lie in so-called "compensatory effects". Compensation can be understood in the sense that more favorable soil conditions (higher soil pH and moisture) outweighed the effects of seemingly stressful abiotic conditions in the form of less available light (higher tree-layer SCA), resulting in high herb-layer species richness. Some of our study sites experienced more favorable light conditions due to the sparser canopy of certain tree species (e.g., Pinus sylvestris) but were at the same time characterized by acidic, nutrient-poor soils. At the other end of the light gradient, our plots included some forest stands dominated by broadleaves or mixed stands with closed tree-layer cover and late-successional tree species able to produce high shade (e.g., Fagus sylvatica, Abies alba). Many of such sites are, however, located on carbonate parent material with base rich calcareous soils, where higher soil pH facilitates plant species richness [114,115]. Shade-tolerant herbaceous perennials that give priority to persistence are the main growth form in the herb layer of these forest ecosystems. In addition, we did not consider plot-level environmental heterogeneity, which might contribute to species coexistence. Spatial variability in resources can be even more important for plant species richness than average site conditions [116,117]. Resource heterogeneity exerts a greater influence on vegetation diversity at the later stages of succession, i.e., in closed canopy mature forest stands [21]. Soil heterogeneity seems to be especially important for typical forest plant species [118]. Moreover, sun flecks can have a prominent effect on fine-scale patterning of herb-layer composition and diversity [119].
Our general observation regarding soil properties was that the organic soil layer was more important for herb-layer vegetation than mineral soil layers. Herb-layer species richness was positively related to increasing soil pH of the organic layer (similar to the work of [22] and [110]), a pattern which confirms that forest sites on limestone and dolomite bedrock with more calcareous soils (characterized by higher pH and cation exchange capacity) host higher plant diversity [114]. The authors of [120] reported that understory plant richness decreased significantly with increasing soil acidity. Another relation to soil properties in our plots was observed for herb-layer species evenness, which was negatively correlated with total nitrogen content in organic soils. This could suggest that herb-layer communities on more productive sites tend to be characterized by the prevalence of competitive nutrient-demanding species with vigorous growth, resulting in uneven distribution of abundances among dominant and subordinate species.
Climatic parameters and topographic and structural factors were, compared to forest stand characteristics and soil properties, relatively unimportant for taxonomy-based diversity measures of herb-layer communities. We assume that the climatic signal was weakly expressed due to the usage of broad-scale and interpolated climatic data rather than in-situ measurements. In the case of herb-layer species richness, the cover of rocks on the surface was positively associated with richness. This can be also interpreted in relation to bedrock and soil conditions as forest sites on carbonate parent material usually show higher surface rockiness. Focusing on herb-layer cover, moss cover on the forest soil was a significant predictor with a hump-shaped trend. This indicates that mosses can be both facilitative (creating favorable microsites for seedling establishment) and inhibitive (competition for growing space and resources) for the cover of herb-layer vascular plants. The only climatic parameter that was included in the final GAM for herb-layer cover was temperature seasonality, but even in this case, the effect was marginally significant.

Trait-Environment Relationships
Trait-based approaches increase universality and allow comparisons between different regions of the globe. Thus, we assume that our results related with functional diversity measures are likely more comparable with other studies [121,122] than findings about taxonomy-based measures.
Analyzing convergence and divergence in plant traits along ecological gradients [47,49] complements patterns observed on the taxonomic (species) level. The understanding of trait variation along environmental gradients is one of the main research directions in plant functional ecology [19]. We included three groups of plant functional traits: clonal traits, bud bank traits and traits belonging to the leaf-height-seed (LHS) scheme [80]. While above-ground LHS traits have been frequently studied, below-ground traits remain largely unexplored, and limited evidence exists in the scientific literature regarding forest herb-layer vegetation [41].
The diversity of clonal traits (clonality, clonal index, type of clonal growth organ, persistence of clonal organ) was most significantly driven by stand characteristics (shrublayer richness), climatic parameters (precipitation seasonality) and soil properties (soil pH of the organic layer). The authors of [41] found that climate had a pervasive effect on determining patterns of clonal trait composition (not diversity) in Italian forest understories. The authors of [49] concluded that clonal traits were significantly more diverse in habitats that represented high soil reaction. Our results agree with this observation as we showed a shift from clonal trait convergence to divergence along the organic soil pH gradient. An interesting pattern was observed for soil texture. On sites with a higher proportion of clay in the upper mineral soil layer, clonal traits tended to converge (i.e., low dissimilarity among species regarding the clonal traits). One possible explanation could be that belowground clonal organs have limited space for growth in soils where clay particles are more abundant, and soils are more compacted.
The diversity of bud bank traits (bud bank size and bud bank depth) was mainly related to altitude, organic soil pH, tree-layer SCA and, similar to clonal traits, to precipitation seasonality. However, bud bank trait diversity changed differently along the precipitation gradient compared to clonal traits. Furthermore, it is worth noting that altitude (here assigned to the group of topographic and structural factors) is actually a variable very closely associated with mean climatic conditions. Across our dataset, altitude exhibited a strong negative correlation with mean annual temperature (r = −0.94) and a positive correlation with precipitation (r = 0.43). Elevation played a major role in shaping plant diversity in the study of [42], where high-elevation plots exhibited reduced functional diversity of below-ground storage organs that may promote cold-tolerance. Our findings similarly suggest that the diversity of bud bank traits declined with altitude, meaning that bud bank trait divergence occurred at low-altitude sites with higher temperatures. This corroborates with the results of [49], who reported that bud bank traits were more diverse in warmer habitats. In addition, our results are also in line with theirs in terms of the dependence of below-ground traits on light conditions. Convergence of bud bank traits occurred in plots with higher tree-layer SCA, meaning that in forest understories with limited light availability, species seem to use similar clonal strategies to cope with shade.
Compared to clonal and bud bank traits, LHS traits were significantly associated with only two soil properties: soil pH and nitrogen content of the organic layer. In their study of LHS trait composition (not diversity) across understory plant communities in Italian forests, the authors of [40] summarized that climate-related factors as well as soil nitrogen and phosphorus availability were the explanatory variables most correlated to trait variation. They explicitly highlighted the need to integrate soil properties as local drivers of above-ground trait diversity. In our case, the diversity of LHS traits decreased along the soil pH and productivity gradient. This agrees with the work of [49] who found a negative correlation between soil reaction and LHS trait diversity over a large set of temperate vegetation types. Furthermore, we have additionally observed an interesting pattern regarding this trait group across our monitoring plots. While clonal and bud bank trait diversity were positively associated with herb-layer species richness, for LHS traits the relationship was negative (results not shown). This suggests that species-poor herblayer communities tend to be functionally more diverse than species-rich communities, pointing to the conclusion that species richness is not always a good surrogate for functional diversity. Moreover, stronger filtering of LHS traits occurred in environments (forest sites) with sufficient nutrient levels but are generally less favorable in terms of light availability, such as dense beech or mixed broadleaf forests on more fertile soil types.

Conclusions
This work investigated the effects of ecological drivers on the vascular plant diversity of herb-layer communities in Slovenian forests. The main strengths of our study are (1) the use of harmonized ICP Forests vegetation records accompanied with detailed in-situ measurements of soil properties and estimations of a wide variety of ecological factors as potential predictors and (2) the inclusion of taxonomy-and trait-based measures of diversity. Moreover, within each diversity facet, three different response variables were distinguished. To date, only a few studies (e.g., [42]) have simultaneously tested the impacts of environmental drivers on both taxonomic and functional aspects. Most studies have focused on comparing species (trait) richness and (functional) composition (e.g., [5,120]) rather than on integrating multiple facets of diversity. Disentangling the relative importance of ecological drivers for different organizational levels can elucidate how abiotic constraints influence plant community assembly and potentially have implications for forest biodiversity conservation and management of forests.
Variance partitioning analyses revealed the prevalent effect of forest stand characteristics and soil properties on the majority of investigated taxonomic and functional trait variables. Some common patterns emerged. One evident example was the lack of a climatic signal for taxonomy-based measures of diversity. Although we can draw such general conclusions based on the comparison between taxonomy-and trait-based diversity measures, it is important to note that even within the same diversity facet (taxonomic and functional), each response variable responded quite differently to ecological factors. Our next research step will be dedicated to the quantification of how species (e.g., using ordination methods) and functional composition (e.g., by combining the fourth-corner and the RLQ method; [121,123]) of herb-layer communities are related to ecological drivers. We encourage researchers from other countries involved in the ICP Forests monitoring network to take a similar multi-faceted approach by systematically including functional diversity instead of taxonomic metrics only and maintain that this should be incorporated into biodiversity monitoring plans [42]. Comprehensive comparisons will not only reveal factors influencing species and functional trait diversity but can also provide useful insights related to relevant drivers of temporal dynamics in European forest plant communities, particularly in the face of global environmental change and the biodiversity crisis [4,24,117].  Institutional Review Board Statement: Not applicable (study does not involve human or animals).

Informed Consent Statement:
Not applicable (study does not involve human).
Data Availability Statement: Datasets generated and/or analyzed during the current study are available from the corresponding author on request.