Changes in Species and Functional Diversity of the Herb Layer of Riparian Forest despite Six Decades of Strict Protection

: The herb layer of temperate forests contributes to long-term forest ecosystem functioning and provisioning of ecosystem services. Therefore, a thorough understanding of its dynamics in the face of environmental changes is essential. This paper focuses on the species and functional diversity of the herb layer of riparian forests to verify how these two community components changed over time and under strict protection. The understory vegetation was surveyed on 42 semi-permanent plots in three time periods between 1960 and 2020. The overall pattern in vegetation changes that related to species richness and diversity, functional structure, and habitat conditions was analyzed using ordination and permutation techniques. We found signiﬁcant changes in species composition and the functional structure of herbaceous vegetation over the last six decades. Forests were enriched with nutrient-demanding and alien species. A signiﬁcant increase in functional diversity and the proportion of species with high SLA and canopy height was also observed, whereas changes in habitat conditions were insigniﬁcant. The observed trends indicate that the strict protection of forest communities within small and isolated reserves does not fully protect their species composition. Forest reserves should be surrounded by unmanaged forests and spatially connected to allow species mobility.


Introduction
In temperate forests, most of their biodiversity is concentrated in the herb layer [1]. Thus, processes taking place in this vegetation stratum may contribute to long-term forest ecosystem functioning and, at the same time, affect the ability of forests to provide ecosystem services [2][3][4]. However, currently, forests face severe threats which can act independently or synergistically and, most importantly, affect forests globally [5]. One of these risk factors is climate change, reflected by rising temperatures [6], extreme weather events such as prolonged droughts [7] or late-spring frost events [8], and intensification of evaporation [9]. These all overlap with rapid thermophilization [5,10], eutrophication [11,12], and plant invasions [13], which change the compositional and functional structure of forests at an unprecedented pace and scale.
Therefore, it is crucial to recognize the current state of forest communities and, more importantly, determine the proportion of biodiversity that can be maintained within forest ecosystems, and the role of protected areas (PAs) in this process. The latter is considered to be one of the primary and effective tools in mitigation of the consequences of rapid changes in a global environment, by reducing both a decline in biodiversity in various types of 1958, we may assume that they have been subjected mainly to the natural processes of forest dynamics. However, due to the observed global trends (invasions of species, warming, eutrophication) and the small area of studied reserves, we hypothesize that herbaceous vegetation's species composition and functional structure have significantly changed over time, which is consistent with most literature sources.

Study Area
We conducted research in three forest nature reserves established in 1958: "Kanigóra" (5.4 ha), "Grodziska Ryczyńskie" (1.82 ha), and "Zwierzyniec" (8.55 ha). They are located within a sizable riparian forest complex with an area of 2004 hectares [52], preserved in the Odra River Valley. The Odra forms the 30th longest and 17th largest river catchment area of Europe, including Russia (Figure 1), and is an ecological corridor important on the continental scale [53]. The entire area is located on the Silesian Lowland, in the southwestern part of Poland. The elevation varies between 120 and 140 m a.s.l., and the plant communities are typical of Central European lowlands.
Forests 2022, 13, x FOR PEER REVIEW 3 of structure over time; (ii) to recognize which environmental factors may be responsible f the observed changes. As the studied forests have been continuously protected as natu reserves since 1958, we may assume that they have been subjected mainly to the natur processes of forest dynamics. However, due to the observed global trends (invasions species, warming, eutrophication) and the small area of studied reserves, we hypothesi that herbaceous vegetation's species composition and functional structure have signi cantly changed over time, which is consistent with most literature sources.

Study Area
We conducted research in three forest nature reserves established in 195 "Kanigóra" (5.4 ha), "Grodziska Ryczyńskie" (1.82 ha), and "Zwierzyniec" (8.55 ha). Th are located within a sizable riparian forest complex with an area of 2004 hectares [52 preserved in the Odra River Valley. The Odra forms the 30th longest and 17th largest riv catchment area of Europe, including Russia (Figure 1), and is an ecological corridor im portant on the continental scale [53]. The entire area is located on the Silesian Lowland, the southwestern part of Poland. The elevation varies between 120 and 140 m a.s.l., an the plant communities are typical of Central European lowlands. The soils in the study area are alluvial and formed during river floods. Gley so dominate in the northern and western parts, while brown soils occur in the central an eastern parts, showing high compactness, acidity, and high moisture [52].
The mean annual temperature is typical of the warmest region in Poland and is abo 8 °C. The warmest month is July (mean 18 °C), while the coldest is January (mean 1.1 °C The total annual rainfall varies between 560 and 580 mm. Snow cover in this region do not last long; its average durability is the lowest in Poland and amounts to less than days a year. The growing season lasts more than 230 days, and it has gradually lengthene as a result of climate change [54].
The studied forests are one of the critical objects of protection in the Natura 2000 si PLH020017 "Grądy w Dolinie Odry" covering the area of 8756.34 ha [55]. Initially, mosaic of riparian, oak-hornbeam and alder forest [56][57][58][59] were recognized her The soils in the study area are alluvial and formed during river floods. Gley soils dominate in the northern and western parts, while brown soils occur in the central and eastern parts, showing high compactness, acidity, and high moisture [52].
The mean annual temperature is typical of the warmest region in Poland and is above 8 • C. The warmest month is July (mean 18 • C), while the coldest is January (mean 1.1 • C). The total annual rainfall varies between 560 and 580 mm. Snow cover in this region does not last long; its average durability is the lowest in Poland and amounts to less than 40 days a year. The growing season lasts more than 230 days, and it has gradually lengthened as a result of climate change [54].
The studied forests are one of the critical objects of protection in the Natura 2000 site PLH020017 "Grądy w Dolinie Odry" covering the area of 8756.34 ha [55]. Initially, a mosaic of riparian, oak-hornbeam and alder forest [56][57][58][59] were recognized here. However, recently they naturally transformed into forests of lower dampness, dominated by hornbeam stands, Stellario holosteae-Carpinetum betuli Oberdorfer 1957 [60]. Outside the area of the reserves, forests are intensively managed. Forest management in the area is carried out by the Oława Forest District, which obtains at least 92,000 m 3 of wood annually, mainly coniferous wood, but also oak, ash, and alder. Tourist and recreational pressures are also intense, especially in the southern part of the area which is crossed by many forest roads.

Sampling
The previous surveys were carried out twice: in 1958-1961 [56] (hereafter 1960) and 1988-1991 (hereafter 1990) [57][58][59]. The most recent survey was conducted in 2019-2020 (hereafter 2020), covering the whole vegetation season, similar to the above-mentioned studies. The location of the most recent plots was determined based on the descriptions of the locations included in the study of Kuczyńska [56], undertaken in 1960, considered as a baseline survey (see Supplementary Materials, Table S1). Therefore, we treated the newly established plots as semi-permanent. In total, we resurveyed 42 semi-permanent plots: 20 in the Zwierzyniec, 7 in the Grodziska Ryczyńskie, and 15 in the Kanigóra reserves. The location of each plot was determined using GPS coordinates (see Supplementary Materials, Table S1). The plots' size was normalized to 100 square meters. According to the Braun-Blanquet approach, vegetation composition was sampled in the form of phytosociological relevés [61], with abundance determined for each species in a separate vegetation layer (tree, shrub, and herb layers) on a seven-point scale. The most recent survey was conducted with permission of the Regional Directorate for Environmental Protection in Wrocław (decision number WPN.6205.58.2019.AR and WPN.6205.39.2020.MR).
The nomenclature of vascular plants follows Euro + Med PlantBase [62]. Mosses were omitted due to random data from the past periods of investigation.

The Overall Pattern in Vegetation
The differences in species composition of the herb layer between surveys were analyzed using JUICE 7.1.21 software [63] with a modified TWINSPAN algorithm [64]. All species occurring in the herb layer, including herbaceous species and saplings of shrubs and trees, were considered in these analyses.
Differentiating species for the analyzed surveys were determined using the Φ [phi] coefficient as a measure of fidelity [65]. Species with Φ ≥ 25 (0.25 × 100), constancy ≥ 20%, constancy ratio [66] higher than 1.5 and significant concentration in a particular cluster, tested by the Fisher's exact test (p < 0.05), were considered to be diagnostic (i.e., reaching the highest fidelity and frequency) for particular survey periods. In addition, the total number of species (species richness) and the Shannon-Wiener Index were calculated using a JUICE function for each relevé in the subsequent survey times. In this respect, differences between the surveys were tested using a Kruskal-Wallis test and Dunn's post hoc test. Moreover, in the JUICE software, we calculated the frequency and coverage of alien species in vegetation plots referring to each studied survey. The analyses were performed jointly on the entire dataset. Alien species were determined following Tokarska-Guzik et al. [67].

Habitat Conditions
Due to the lack of direct measurements of environmental conditions in the baseline and current research, we estimated habitat conditions indirectly using litter quality index (LQ) and Ellenberg's Indicator Values (EIVs). Both measures derived from the species composition of studied plots.
The change in the composition of the canopy, including tree and shrub layers, was quantified using mean LQ. The LQ ranges from 1 (slow litter decomposition rate) to 5 (fast litter decomposition rate). Indices scores for tree and shrub species follow those used by Verheyen et al. [68] and Šipoš et al. [69] (see Supplementary Materials, Table S2). We calculated the mean value of LQ weighted by species cover for shrub (LQs) and tree (LQt) layers separately for each plot. Differences between the LQ values between subsequent surveys were tested using the Kruskal-Wallis test and Dunn's post hoc test. The EIVs, weighted by species cover, were calculated for light, temperature, moisture, continentality, soil reaction, and nutrients. We used Ellenberg's Indicator table [70], modified and implemented in JUICE software [71].
We conducted a detrended correspondence analysis (DCA) to explore differentiation patterns in plant species composition of plots between the subsequent surveys and to check the percentage of variation explained. For this purpose, species cover was transformed with Log (1 × X + 1) and rare species were down-weighted. The distribution of the sample groups on the DCA diagram was visualized and interpreted by the passive plotting of the Year supplementary variable. The DCA was conducted in CANOCO 5.12 software (Microcomputer Power: Ithaca, NY, USA) [72].
To identify the statistical significance of correlations between the DCA sample scores and mean randomized EIVs (treated as proxy data to environmental changes) for relevés, a modified permutation test with 499 unrestricted permutations was conducted. The test was performed with MoPeT_v1.2.r script [73] in R software [74]. The Spearman rho coefficient expressed the correlation between EIVs and sample scores. Permutational analysis of variance (one-way ANOVA on the mean randomized EIVs) and modified permutation test (with 499 unrestricted permutations) were also calculated using MoPeT_v1.2.r to determine which EIVs differentiate the subsequent surveys. Using permutational ANOVA is an alternative to other tests under non-normal conditions, because it does not operate under the assumption of normality and uses actual scores [75]. This analysis allows determining possible changes in the environmental conditions influencing the forest floor over time.

Plant Functional Traits
In this study we focus on the plant functional traits which are considered the core of plant life strategy [76] and belong to response traits [77]. These are specific leaf area (SLA, mm 2 ·mg −1 ; the ratio of leaf area to leaf dry mass), leaf area (LA, mm 2 ; one-sided or projected area of an individual leaf), plant height (H, m; plant vegetative height at maturity), and seed mass (SM, g).
Higher values of the SLA reflect higher investment in photosynthesis and fast growth. They are typical of short-lived species which do not invest in support and storage organs, and often are weakly protected against herbivores [2]. At the same time, their response to the spatial patchiness of light and soil resources may be more flexible [76]. On the other hand, most of the species with low SLA are long-lived and marked by slow growth, high investment in support and storage organs, and higher resistance to drought and herbivory [2]. Therefore, an increase in SLA values may suggest a higher level of competition, but simultaneously a lower level of stressful conditions [42,78]. Species with higher LA are expected to dominate in shadier, more humid habitats. However, any stressful condition such as heat stress, cold stress, drought stress, or nutrient stress tends to select for relatively small leaves [79]. Plant height is associated with a competitive ability and reflects trades off between avoidance and tolerance of environmental conditions (e.g., nutrient and light availability) [42,76]. Seed mass expresses the chance of species for successful dispersion and establishment in a new area [76]. Additionally, it provides information on seed bank longevity which is expected to be lower in habitats with a high level of environmental stress [80]. It has already been proved that the SM increases with the shadiness of the habitat and mean vegetation height. SM is also expected to decrease in more disturbed habitats [81].
We obtained trait values from the LEDA database [82]. We used mean species values for traits and did not consider intraspecific trait variation [43,83], (see Supplementary Materials, Table S3). We excluded all woody species (tree, shrub) from the dataset in this part of the data analysis. Therefore, only herbaceous plants (listed in Supplementary Materials, Table S4) were assigned mean trait values. Trait attributes were analyzed for 92 out of the 97 species recorded in the herb layer, whereas the remaining five species, without full representativeness of trait attributes, were removed from the dataset [84]. The functional structure of the plant community can be reflected by different metrics [85][86][87]. In this study we focused on the community weighted mean of traits (CWM) [85,88] and Rao's quadratic entropy (FD Q ) [89]. The former is widely used in vegetation ecology and is a mean of trait values present in the community weighted by the relative abundance of taxa bearing each value [85]. In other words, the CWM index reflects the trait value of the dominant species in a community [42]. The FD Q is an index of functional diversity which incorporates species abundances and a measure of the pairwise functional differences between species [89]. It quantifies how much dominant species diverge in their trait values; thus, its increased value may indicate a higher degree of niche differentiation [41,42]. Moreover, it can be calculated for more than one trait and is not affected by species richness [89]. The latter feature of this index is critical because, in this study, we analyzed vegetation within semi-permanent plots surveyed by different researchers.
Differences in CWMs and FD Q between subsequent surveys were tested using the Kruskal-Wallis test and Dunn's post hoc test. Additionally, in CANOCO 5.12 software [72] we conducted a redundancy analysis (RDA) in which the explanatory variable was the year of each survey and CWMs were treated as response variables. All traits we used were continuous in character and log-transformed before the analyses.
out of the 97 species recorded in the herb layer, whereas the remaining five species, without full representativeness of trait attributes, were removed from the dataset [84].
The functional structure of the plant community can be reflected by different metrics [85][86][87]. In this study we focused on the community weighted mean of traits (CWM) [85,88] and Rao's quadratic entropy (FDQ) [89]. The former is widely used in vegetation ecology and is a mean of trait values present in the community weighted by the relative abundance of taxa bearing each value [85]. In other words, the CWM index reflects the trait value of the dominant species in a community [42]. The FDQ is an index of functional diversity which incorporates species abundances and a measure of the pairwise functional differences between species [89]. It quantifies how much dominant species diverge in their trait values; thus, its increased value may indicate a higher degree of niche differentiation [41,42]. Moreover, it can be calculated for more than one trait and is not affected by species richness [89]. The latter feature of this index is critical because, in this study, we analyzed vegetation within semi-permanent plots surveyed by different researchers.
Differences in CWMs and FDQ between subsequent surveys were tested using the Kruskal-Wallis test and Dunn's post hoc test. Additionally, in CANOCO 5.12 software [72] we conducted a redundancy analysis (RDA) in which the explanatory variable was the year of each survey and CWMs were treated as response variables. All traits we used were continuous in character and log-transformed before the analyses.

Changes in Species Composition
The analysis conducted in JUICE revealed differences in species composition between studied surveys. We identified differentiating species for each survey. The first and the second surveys were marked by species typical of riparian forests. Although these taxa were also found during the recent survey, their frequency and cover in particular plots were significantly lower (Table 1). A distinctive feature of the species composition in 2020 was the presence of the nutrient-demanding and fast-colonizing species, which currently play a diagnostic role within studied forest communities (e.g., Impatiens parviflora, Chelidonium majus, Galeopsis tetrahit, Chaerophyllum aromaticum). The latter period was also marked by the increased proportion of juveniles of Acer campestre, A. pseudoplatanus, Carpinus betulus, Crataegus monogyna, Fraxinus excelsior, Quercus robur, Sambucus nigra, and Tilia cordata. Changes in the proportion of alien species were also clearly visible (Table 2). During the baseline survey (1960) any alien species were recorded. Since 1990 Impatiens parviflora has been observed. However, its frequency and cover (mean cover 2%) were low. In 2020, Impatiens parviflora became one of the dominant species (mean cover 14%), and six other alien species were recorded. All alien species, except for Impatiens parviflora, occur with low frequency and cover (less than 2%). However, some of them, for example Solidago gigantea, are known as highly invasive and may affect the species composition of studied forests in near future.

Temporal Changes in Herbaceous Vegetation Pattern
DCA ordination diagram ( Figure 3) shows a species compositional pattern for all subsequent surveys, with a passive projection of Year variable. The DCA results revealed that the first and the second axes explained 9.64% and 6.81% of the compositional variability of the studied forests, respectively.

Temporal Changes in Herbaceous Vegetation Pattern
DCA ordination diagram ( Figure 3) shows a species compositional pattern for a subsequent surveys, with a passive projection of Year variable. The DCA results revealed that the first and the second axes explained 9.64% and 6.81% of the compositiona variability of the studied forests, respectively. The results of DCA suggest that the distribution of the sampled plots correspond with the vector of time, indicative by a shift from the lower left (1960) towards the uppe right (2020) of the ordination diagram ( Figure 3). The same diagram with species instead of sample plots showed that the most recent survey is marked by the nutrient-demandin species, such as Impatiens parviflora, Stellaria media, Galium aparine, and additionally som juveniles (Figure 4). The results of DCA suggest that the distribution of the sampled plots correspond with the vector of time, indicative by a shift from the lower left (1960) towards the upper right (2020) of the ordination diagram ( Figure 3). The same diagram with species instead of sample plots showed that the most recent survey is marked by the nutrient-demanding species, such as Impatiens parviflora, Stellaria media, Galium aparine, and additionally some juveniles (Figure 4).
Correlation analysis of mean EIVs and main compositional gradients (reflected by the DCA sample scores) showed statistically significant correlations only along the second DCA axis. The distribution of sample scores was correlated to the EIVs for moisture (p < 0.02) and nutrients (p < 0.05) ( Table 3).  Correlation analysis of mean EIVs and main compositional gradients (reflected by the DCA sample scores) showed statistically significant correlations only along the second DCA axis. The distribution of sample scores was correlated to the EIVs for moisture (p < 0.02) and nutrients (p < 0.05) ( Table 3).   On the other hand, ANOVA revealed no statistically significant differences in the mean EIVs between the studied time periods (Figure 5), although the EIV for light showed a steady increase over the last 60 years. On the other hand, ANOVA revealed no statistically significant differences in the mean EIVs between the studied time periods (Figure 5), although the EIV for light showed a steady increase over the last 60 years.

Changes in Litter Quality Index
We found no statistically significant changes in the litter quality index for trees (LQt) between 1960 and 2020 (p = 0.25), although there was a difference between the first and the last surveys and the middle one (in both cases p < 0.001, see Figure 6a). On the other hand, the litter quality index for shrubs (LQs) has gradually decreased (see Figure 6b), with statistically important differences between 1960 and 2020 (H = 11.4, p < 0.002).

Changes in Litter Quality Index
We found no statistically significant changes in the litter quality index for trees (LQt) between 1960 and 2020 (p = 0.25), although there was a difference between the first and the last surveys and the middle one (in both cases p < 0.001, see Figure 6a). On the other hand, the litter quality index for shrubs (LQs) has gradually decreased (see Figure 6b), with statistically important differences between 1960 and 2020 (H = 11.4, p < 0.002).

Changes in Functional Structure
Analysis of the CWMs of functional traits suggests that some have been unchanged in the herb layer since 1960. For two functional traits, we found changes with a high level of significance between 1960 and 2020. Changes in CWM for SLA were non-significant

Changes in Functional Structure
Analysis of the CWMs of functional traits suggests that some have been unchanged in the herb layer since 1960. For two functional traits, we found changes with a high level of significance between 1960 and 2020. Changes in CWM for SLA were non-significant between 1960 and 1990 (Figure 7a), but during the last 30 years the value of this trait suddenly increased, with a high level of significance (H = 54.76, p << 0.001). On the other hand, the mean height of the herb layer significantly increased between 1960 and 1990 (p << 0.001) (Figure 7b), and between 1960 and 2020 (p < 0.001). Figure 6. Summary box-and-whisker plots with jitters presenting changes in litter quality of trees (LQt) (a) and shrubs (LQs) (b) between three studied surveys of the riparian forest in the Odra River Valley (SW Poland). The central line of each box indicates the median value, box boundaries the lower (25%) and upper (75%) quartiles, and whiskers the range of values.

Changes in Functional Structure
Analysis of the CWMs of functional traits suggests that some have been unchanged in the herb layer since 1960. For two functional traits, we found changes with a high level of significance between 1960 and 2020. Changes in CWM for SLA were non-significant between 1960 and 1990 (Figure 7a), but during the last 30 years the value of this trait suddenly increased, with a high level of significance (H = 54.76, p << 0.001). On the other hand, the mean height of the herb layer significantly increased between 1960 and 1990 (p << 0.001) (Figure 7b), and between 1960 and 2020 (p < 0.001).  There were no significant differences between the first and the last period of investigation concerning LA (p = 0.13, Figure 7c) and SM (p = 0.4, Figure 7d).
The results of the RDA also suggested that the CWMs of studied plant functional traits were not equally responsible for the observed changes in functional structure of individual patches. All plots investigated between 1960 and 1990 showed a high consistency of functional features. Only some of them were marked by a high level of CWMs for seed mass and height. However, in the most recent survey, several analyzed patches showed much higher values of SLA, and partially also of height, which has been already confirmed by the presented analysis of changes in functional traits over time (Figure 8). The first axis, that explained 17.4% of the variation, was correlated with changes in SLA over time. The second axis, explaining 38.38% of the variability, seemed to correlate with leaf area, height, and seed mass, accounting more for the internal variability of the investigated plots ( Figure 8). In the first period of investigation plots were marked by the higher values of the CWMs for LA and SM, and during the second and third surveys the CWMs for H increased. already confirmed by the presented analysis of changes in functional traits over time (Figure 8). The first axis, that explained 17.4% of the variation, was correlated with changes in SLA over time. The second axis, explaining 38.38% of the variability, seemed to correlate with leaf area, height, and seed mass, accounting more for the internal variability of the investigated plots (Figure 8). In the first period of investigation plots were marked by the higher values of the CWMs for LA and SM, and during the second and third surveys the CWMs for H increased. Functional diversity FDQ (measured as Rao's quadratic entropy) of the investigated forests were stable from 1960 to 1990 (p = 0.35), and after 1990 it rapidly increased ( Figure  9). There were significant differences between 1960 and 2020, as well as between 1990 and 2020 (p < 0.001). Functional diversity FD Q (measured as Rao's quadratic entropy) of the investigated forests were stable from 1960 to 1990 (p = 0.35), and after 1990 it rapidly increased ( Figure 9). There were significant differences between 1960 and 2020, as well as between 1990 and 2020 (p < 0.001).

Species Composition and Diversity Change
We observed significant qualitative and quantitative changes in the herbaceous

Species Composition and Diversity Change
We observed significant qualitative and quantitative changes in the herbaceous vegetation of the studied riparian forests over the last six decades, which refer to species composition, richness, and diversity. Our results revealed that the species composition of the studied forests has shifted from the community with a significant proportion of forest specialists [90] (e.g., Aegopodium podagraria, Ajuga reptans, Paris quadrifolia, Symphytum tuberosum) to the one in which nutrient-demanding species play a differentiating role (Impatiens parviflora, Chelidonium majus, Galeopsis tetrahit) ( Table 1).
It is also worth mentioning that in these forests there are neophytes which considerably increased their frequency and coverage over the study time ( Table 2). Their presence agrees with earlier studies which indicated that the riparian forests are the communities highly vulnerable to plant invasions [29,91,92]. The observed compositional change translates into alterations in species richness and diversity-both metrics showed a significant increase in their values, especially during the last three decades (Figure 2). Therefore, studied forests follow the widely observed pattern of changes currently taking place in forest ecosystems of the temperate zone. The replacement of species typical of forests by disturbance-adapted and nutrient-demanding taxa has been recently extensively reported [12,68,93].
Interestingly, this phenomenon has not yet caused a decline in species richness at the local scale. The loss or decline in the proportion of forest specialists is balanced by species gains [93], which is also true for the present study. The gradual increase in species richness that we observed may be explained by the intermediate disturbance hypothesis, described by Connell [94], who stated that disturbances at an intermediate level maximized species richness. However, an increase in this parameter is often due to the appearance of early-successional, light-demanding, and invasive species that often replace shade-tolerant forest specialists [95]. For these reasons, the total species number per se cannot be treated as the only indicator of forest conservation values, since such an increase is often caused by the entry into the forest by species responding to local disturbance [96]. Our study shows the same patterns, since these gains in species number are associated with the increased role of species that, over a long period, may negatively affect the pool of true forest species associated with forest interiors.

Possible Reasons for the Change in Species Composition
According to Verheyen et al. [68], the increasing proportion of nutrient-demanding species may result from a rising density of tree canopy and changes in its composition in favor of tree species with a higher LQ index [68]. However, in the present study, LQ remained similar over time (LQt), or decreased (LQs); thus, we may assume it did not affect the species turnover. The increase in species richness from the initial study up to the recent survey may also be connected with high levels of nutrients, which is in line with findings of Bernhardt-Römermann et al. [97] and Naaf and Kolk [98], but in contrast to Verstraeten et al. [99] who revealed a decline in the mean species richness on neutral-sites. Regarding the presence of alien species in PAs, studies conducted in Central Europe pointed out that nature reserves were also affected by the invasion of species, although, in general, it seemed to be less harmful than in the case of managed forests [26,29,100]. Notably, the study of Pyšek et al. [34] confirmed that the older the nature reserve is, the more resistant to invasions. Such tendencies can be due to the usually higher species richness observed in nature reserves, and the much better condition of the plant communities protected inside them [100]. It was also examined that, in more diverse and well-preserved (undisturbed) forest communities, alien species coexist with the native ones rather than compete, contributing to higher species richness and diversity, or the lack of significant, qualitative, and quantitative changes in this respect [34,68,101].

Changes in Vegetation Type and Habitat Conditions
Our findings also suggest that although the species composition of the studied riparian forests was less diverse (lower species richness and diversity) in the 1960s, the older plots themselves differed more in terms of habitat conditions. It is visible on the DCA diagram, where the distribution of the studied plots along the second axis is correlated with the EIVs for moisture (Figure 3), and nutrients (Table 3), and in the DCA plot presenting the distribution of species (Figure 4). Such an arrangement of the sample scores and species corresponds with the tendencies observed in the nature reserves, reported in the earlier studies [56][57][58][59], that emphasized the presence of mainly two types of forest communities different in reference to moisture conditions, i.e., hardwood floodplain forest and oak-hornbeam forest. Nowadays, the forest communities of the studied nature reserves are more abundant in species, and simultaneously, their species composition is more diverse (Figure 2). However, they seem to be more homogenized in terms of these two metrics, which may be connected with the homogenization of edaphic conditions (e.g., ameliorations, the lack of periodic floods, thus slightly lower moisture). The direction of these changes is not surprising since the transition of the European floodplain forests towards drier communities, such as oak-hornbeam forests, has been widely reported [91,92,102]. Nevertheless, the ANOVA conducted for the entire dataset did not reveal any statistically important differences between surveys in the mean values of the EIVs, although some slight upward and downward trends are visible in reference to light and moisture, respectively ( Figure 5). Such results may suggest that studied communities have not substantially changed in terms of habitat conditions over the last 60 years. However, this would be an exception in the background of recent studies, in which temporal changes in environmental conditions are widely reported from deciduous forests [97][98][99]103].

Functional Structure
Studied riparian forests have changed significantly over time in terms of functional structure. The analysis of the CWMs revealed that herbaceous vegetation has become dominated by species with higher values of the SLA and height (H) (Figure 7). Such species, which are usually shade-tolerant, highly competitive, and less stress-resistant (Tables 1 and 2), gradually replace typical forest species with lower values of the abovementioned traits, thus, to a greater extent shaping the present functional structure of the studied forest communities [93].
We also found significant changes in the functional diversity (FD Q ) [89], which has risen over time. Not only was the herb layer of the riparian forests poorer in species, but also the dominant species were more similar to each other regarding their plant traits. The lower values of FD Q may suggest that herbaceous vegetation was subjected to more stressful conditions in the past (e.g., more frequent floods and lower light availability) and species had to compete for resources [41]. Nowadays, studied forests are richer in species that simultaneously differ from each other in life strategies. This may indicate a higher level of niche differentiation and generally more optimal and less stressful conditions for all species in the herb layer [101].
Moreover, Vanneste et al. [4] revealed that older forests might support a higher variation in functional traits (e.g., SLA, H, seed mass) because they provided a wider variety of microhabitats, due to higher abiotic and biotic heterogeneity. Therefore, we can assume that the rise in FD Q presented in this study may result from strict protection, which has enhanced the formation of numerous microhabitats. Thus currently, species in the herb layer coexist rather than compete for resources. However, paradoxically, this seemingly positive change may facilitate the establishment of alien species which avoid competition and may become more abundant in these newly created niches [104].

Conclusions
Our results revealed visible shifts in herbaceous vegetation of studied forest communities during the last six decades. These changes refer not only to species composition but also to the functional structure of surveyed riparian forests. Differences between surveys mainly relate to the increase in the proportion and abundance of nutrient-demanding and alien species, whose presence is a distinctive feature of the most recent plots collected in 2020. We also noted a statistically significant increase in the proportion of species with high SLA and height-features typical of species with a ruderal life strategy [105]. Surprisingly, we did not detect any significant changes in habitat conditions expressed by the EIVs. However, a visible species turnover is reflected by a gradual replacement of forest species with ruderal ones, regardless of that.
On the one hand, the observed trends may suggest that the strict protection of forest communities, within small and isolated reserves, does not fully protect the species combination typical of these communities, against the processes of ruderalization understood not only as changes in the species richness or diversity, but also in the functional structure of the community. This happens even if PAs are located within a large (over 2000 ha) forest complex, but subjected to intense economic (forest management) and tourist pressures. On the other hand, strict protection may have contributed to creating new microhabitats, thus niche differentiation and the reduction of the level of habitat disturbances. These positive changes are reflected by the rise in functional diversity, suggesting that present communities harbor more species that differ in terms of functional traits than the same communities in the past. The problem is that the beneficiaries of these changes are nutrient-demanding, fast-colonizing, and alien species instead of the native ones typical of forest communities.
Therefore, this study shows that a broader perspective is highly needed in the planning and management of the strictly protected areas. It seems that the interior of the protected areas and their surroundings are essential to protect their biotic diversity effectively. This may indicate the necessity to create wide buffer zones around such forest nature reserves, that will be excluded from forest management and whose role will be only to reduce the negative impact of the global and local environmental factors. This approach has been proposed for the first time by Pyšeket al. [34] and is called the Several Small Inside Single Large (SSISL) model, and, in the context of the results presented, it seems to be an appropriate solution to maximize the species richness of nature reserves and to minimize the risk of plant invasions. Therefore, to fulfill the current assumptions of the EU biodiversity strategy concerning the establishment of areas of strict protection, several factors should be taken into account. Firstly, even when protected, forested areas are not extensive in their area (preferably, the bigger they are, and the more diverse in habitat, the better), and they should not be surrounded by intensively managed forests. Secondly, forest nature reserves should be spatially connected to allow free species exchange. Otherwise, when small and isolated, they are prone to undesirable changes in their species composition and functional structure, as we observed.