Field surveys were conducted in 2014, 2015, 2021, 2022, and 2023 across a total of 261 stream sites (Figure 1). A subset of 137 sites was sampled at least twice during the monitoring period (Supplementary Appendix A in Supplementary Material). Macrophytes were sampled from both banks and the channel, when feasible, by wading upstream along a 100-m section of the river reach. Sampling was conducted consistently during the macrophyte growing season (late spring to late summer) across all years, thereby minimizing the influence of seasonal variability on the analyses. Species were identified primarily in the field. Βryophytes and charophytes specimens were collected and transferred to the laboratory for further identification. Abundance data were collected using the qualitative five-point ordinal scale defined in the IBMR methodology (Haury et al., 2006), which is appropriate for calculating IBMR and associated ecological quality classes (Aguiar et al., 2014). However, due to lack of quantitative abundance estimates (e.g. percent cover or individual counts), biodiversity analyses were based on presence–absence data to ensure methodological robustness. Thus, both analyses used the same macrophyte survey dataset, differing only in data treatment (ordinal abundance for IBMR_GR and presence-absence for biodiversity metrics). Furthermore, species were classified according to aquaticity level and floristic group as described in Aguiar et al. (2014). Aquaticity is a qualitative measure of a species’ affinity to water and ranges from 1 to 8. In this study, only species with aquaticity values from 1 (exclusively aquatic plants) to 5 (hygrophilous taxa that may occur submerged during part of the year) were considered. A detailed description of the aquaticity levels is provided in the Supplementary Material (Supplementary Appendix B). Floristic groups were defined as follows: ALG: algae; PTE: pteridophytes; BRM: bryophytes (mosses); BRL: bryophytes (liverworts); PHY: hydrophytic phanerogams; PHE: helophytic phanerogams; PHG: hygrophilous phanerogams; PHX: other phanerogams.

Ecological quality based on macrophyte data was assessed using the IBMR_GR index. The IBMR_GR has been calibrated through the MedGIG intercalibration exercise (Aguiar et al., 2014) and is the official national method for the ecological quality assessment of rivers in Greece based on aquatic macrophytes (Papastergiadou, 2015; 2023). The index was calculated for all sites according to: IBMR = ∑ i E i K i C S i ∑ i E i K i where E i is the coefficient of ecological amplitude for a given species i, K i is the cover scale, and C S i is the species-specific score indicating tolerance to organic pollution. The index is normalized between 0 and 1, and classification into quality classes follows the methodology described in detail by Stefanidis et al. (2022). We then attempted to summarize the transitions in ecological quality for sites that were assessed in multiple years. A dataset of consecutive year pairs was created for each site, recording the change in ecological quality scores. Each record indicates whether the quality improved, declined, or remained stable between the 2 years.

Temporal beta diversity analyses were based on the 137 sites that were sampled at least twice during the study period. To quantify temporal changes in macrophyte community composition, we calculated within site beta diversity using the Jaccard dissimilarity index. Presence–absence data were compiled for each sampling site and year. For every site sampled in at least 2 years, Jaccard dissimilarity was computed between all possible pairs of years. Mean beta diversity values were then summarized for each pair of years to assess overall temporal turnover patterns in the dataset.

In addition, we assessed whether the observed beta diversity of macrophyte communities in each year differed from random expectations. We calculated the standardized effect size (SES) of total beta diversity using a null-model approach based on the Jaccard dissimilarity index. Presence–absence data for all species were organized in a site-by-species matrix for each year. The observed mean Jaccard dissimilarity among all site pairs within a year was computed as a measure of total beta diversity. To generate the null distribution, species occurrences were randomized 999 times by shuffling presences among sites while maintaining the overall frequency of each species. For each permutation, the mean Jaccard dissimilarity was recalculated, producing a null distribution of beta values. The SES was then obtained as: S E S = β o b s − β n u l l ¯ S D β n u l l where positive SES values indicate greater compositional differentiation than expected by chance (community differentiation), and negative values indicate lower differentiation (community homogenization). Values beyond ±1.96 were considered statistically significant. This analysis allowed us to evaluate whether the spatial structuring of macrophyte assemblages within each year reflected deterministic environmental processes or random species turnover (Leprieur et al., 2012; Ding et al., 2024).

To test whether community composition differed significantly among years, we performed a Permutational Multivariate Analysis of Variance (PERMANOVA) using the function adonis () from the vegan R package. The analysis evaluated differences between years as a categorical factor. To identify which specific year pairs differed significantly, we conducted pairwise PERMANOVA (pairwise adonis) tests.

Finally, we evaluated the relative contribution of each sampling unit to overall beta diversity, by calculating Local Contributions to Beta Diversity (LCBD) following Legendre and De Cáceres (2013). LCBD can be used as a metric of ecological uniqueness, highlighting sites with a distinct composition and thereby they have higher contribution to beta diversity. To this end, calculating LCBD may provide insights to help prevent community homogenization (Dubois et al., 2020). The significance of LCBD values was assessed using permutation tests, in which species occurrences were randomly reshuffled among sites to generate a null distribution.

A total of 109 species were recorded across the 137 sampled sites (Supplementary Appendix C in Supplementary Material). Almost half of the total species (n = 51) showed a strong affinity (aquaticity 1–3) for aquatic habitats (Figure 2). The species pool included 29 true hydrophytes and 38 helophytes, along with five aquatic ferns, three macroalgae, and nine bryophyte species. The remaining 25 species are characterized as hygrophytes, plants that typically occur in moist environments, including riverbanks and floodplains.

Average species richness per site did not differ significantly among monitoring years (Figure 3A). Richness was highest in 2014 and 2015 (mean values of 12.53 and 13.61, respectively), while the lowest average richness was observed in 2021 (10.57).

Our analysis provides an overview of the ecological quality assessed using macrophytes from 2014 to 2023. Although not the same sites were assessed each year, the results showed only slight variations in the frequency of Good and High quality classes. The proportion of sites that reached the WFD target ranged from 37% in 2014 and 2023 to 45% in 2021 and 2022 (Figure 4). Sites classified as Moderate quality exhibited larger fluctuations, ranging from 14% in 2014 to 36% in 2023. Differences in the distribution of sites among quality classes across years were tested using a chi-square test of independence. No statistically significant differences were detected (χ² = 23.7, df = 16, p = 0.096). Similarly, the temporal variation of the IBMR_GR values per class category did not show significant differences among the sampling years, indicating an overall stability in ecological quality as assessed using macrophytes (Figure 5).

We also summarised transitions in ecological quality for sites assessed in multiple years by recording changes in ecological quality class between sampling years. We found that 56% of year pairs showed no change, indicating that ecological quality remained stable (Figure 6A). Seventeen percent of the year pairs involved a decline in quality, whereas 37% involved an improvement. Remarkably, 63% of the year pairs between 2015 and 2021 showed an improvement in ecological quality, suggesting an overall enhancement of macrophyte ecological status between the two monitoring periods (Figure 6B). However, between 2021 and 2022, almost 14% of the year pairs exhibited a decline in ecological quality. Additionally, sites sampled in 2023 were not assessed in 2022, and nearly 60% of them were characterized by ecological quality below “Good,” indicating a shift toward poorer quality conditions. Although we did not observe large variations in total quality among the years, there were notable transitions between consecutive years that may reflect shifts in macrophyte communities in response to abrupt changes in anthropogenic stressors.

Analysis of beta diversity comparisons among multiple years indicated low dissimilarity in the consecutive comparisons, with mean β = 0.29 for 2014–2015 and 0.29 for 2021–2022, implying relatively stable macrophyte assemblages across adjacent years. In contrast, comparisons among larger temporal gaps showed much higher turnover (Figure 3B). Mean Jaccard dissimilarity was 0.64 between 2014 and 2021, 0.74 between 2015 and 2022, and 0.75 between 2014 and 2022. These results indicate substantial compositional change between the early (2014–2015) and later (2021–2023) sampling periods. Between 2021 and 2023 (2 year gap) Jaccard dissimilarity was 0.42, an intermediate value, suggesting some ongoing change after 2021 but less pronounced than the longer gaps. These results are confirmed by the pairwise PERMANOVA which showed significant temporal differences in macrophyte composition for all year comparisons, except for 2014–2015 and 2021–2022. Overall, this pattern implies temporal stability at short lags and a substantial change in community composition over longer intervals, which could reflect long term ecological shifts due to cumulative environmental pressures or management changes.

The mean values of Local Contribution to Beta Diversity (LCBD) showed clear temporal variation across the study years (Figure 7A). In 2014, the average LCBD was relatively high, followed by a noticeable decrease in 2015. Remarkably, LCBD values were lowest in 2021 and increased slightly in 2022. The highest mean LCBD was recorded again in 2023. These results indicate that the ecological uniqueness, measured as contribution of local sites to beta diversity, fluctuated considerably over time, potentially reflecting changes in environmental or ecological conditions influencing community structure.

In addition, the standardized effect size (SES) of total beta diversity, calculated using Jaccard dissimilarities, varied among sampling years, also revealing temporal shifts in community structure. In 2014, 2015, and 2021, SES values were significantly negative (−2.06, −3.03, and −2.11, respectively), indicating that macrophyte communities were more compositionally similar than expected by chance, consistent with biotic homogenization (Figure 7B). The 2022 assemblages did not differ significantly from null expectations (SES = −0.57), suggesting no strong deviation from stochastic patterns. In contrast, the 2023 communities exhibited a significantly positive SES (+2.26), indicating that sites were more dissimilar than expected under random assembly, reflecting biotic differentiation. Overall, these results reveal a shift from homogenization in earlier years to marked differentiation in 2023, signaling a clear reversal in community assembly patterns.

All EU Member States have designed and implemented mitigation measures targeting key drivers of environmental degradation, primarily pollution and hydromorphological modifications, since the adoption of the Water Framework Directive (2000/60/EC). However, by 2015, only 53% of water bodies across the EU had reached the target of “good” status, and evidence from the second round of river basin management plans suggests that there has been limited improvement in ecological quality in recent years (Giakoumis and Voulvoulis, 2018). Reports from the EEA confirm that the percentage of surface water bodies achieving a “good” ecological status has remained largely unchanged over the past two decades (European Commission, 2025).

The situation in Greece is similar, although the country is characterized by a higher proportion of water bodies with Good or High status compared to other regions of Europe (European Environment Agency, 2025). Nonetheless, no substantial improvement has been observed over time, and there is evidence suggesting that the implemented measures are not fully effective (Figure 4). For instance, a study by Stefanidis et al. (2021b) showed that the probability of a water body within a Natura 2000 site achieving “good” status is similar to that of a water body outside such a site, raising concerns about the effectiveness of protection and management measures within protected areas. In Greece, the implementation of the WFD began with significant delays (Skoulikidis et al., 2021) with consequent difficulties in monitoring and assessment. The first monitoring results cover the period 2012–2015, while the second phase of the monitoring program continued from 2018 to 2023, following a 2 year delay. Despite the caveats in implementation of a systematic monitoring, the results from these two periods indicate only limited improvement over time.

In this study, macrophyte-based assessments confirm this trend and demonstrate that results from macrophytes are consistent with those obtained using other Biological Quality Elements (BQEs). Overall, these findings reinforce the consensus that improvements in ecological status are problematic and that further actions are required to reverse this situation. The causes of stagnation in ecological status are not fully understood, but we can identify three main possible explanations: (i) inadequate implementation of the required measures (ii) climate change (Free et al., 2024), which shifts management targets and hinders the effectiveness of measures, and (iii) co-acting stressors that, through synergistic or antagonistic effects, may cause unexpected ecological changes in water bodies that cannot be anticipated by managers (Lemm et al., 2021). Obviously, further research is required to precisely locate the reasons for the lack of quality improvement in Greece, although it is likely that a combination of all three aforementioned factors is present in many regions.

Our results showed only small changes in ecological quality as assessed with macrophytes (Figure 4). Normally, this finding would suggest that communities have remained largely unchanged, exhibiting limited differentiation over time that reflects an unaltered ecological quality. Here, however, this was only partially true. Although no significant variation was observed among consecutive monitoring years, noticeable changes were evident over the decade. Beta diversity among consecutive years was markedly lower than beta diversity among larger time gaps. The standardized effect size (SES) of total β-diversity varied markedly among years (Figure 7). During 2014–2021, SES values were significantly negative, indicating community homogenization relative to null expectations. In contrast, 2023 exhibited a significantly positive SES, suggesting strong differentiation in macrophyte composition. Consistent with this, a shift in LCBD values was also apparent in 2021–2022 with higher values recorded in 2023 indicating an increase in the ecological uniqueness and reduced biotic homogenization (Figure 7). Overall, these findings suggest that 2023 represents a coordinated shift in community composition caused by either environmental change or anthropogenic interventions (e.g. management measures) rather than random variability. What is also important to note is that average species richness remained practically unchanged throughout the monitoring period whereas beta diversity increased (Figure 3). Studies that examine temporal dynamics of alpha and beta diversity in freshwater ecosystems have also reported various patterns of change with alpha diversity either decreasing or increasing and beta diversity increasing (Mi-Jung and Eui-Jin, 2024; Liao and Soininen, 2025). These findings underscore the importance of including beta diversity assessments in biodiversity research, as alpha diversity on its own is not sufficient to capture biodiversity responses to ecological change.

Furthermore, temporal beta diversity can reveal patterns of community homogenization, indicating when communities are becoming more similar. Investigating biotic homogenization is challenging and requires extensive data and resources, but several studies so far have identified cases of reduced heterogeneity among biotic communities in freshwater ecosystems. For example, in the Lower Danube Delta, a study spanning nearly 20 years of macrophyte monitoring revealed significant signs of community homogenization accompanied by species richness loss (Beracko et al., 2025). In lakes of Minnesota, homogenization of macrophyte communities was shown to be driven by invasive species, highlighting the importance of regional mitigation measures (Muthukrishnan and Larkin, 2020). However, homogenization appears to be region-specific, as other long-term studies spanning several decades have not reported any loss of biotic uniqueness or spatial β-diversity. For instance, in boreal lakes, Lindholm et al. (2019) found no evidence of macrophyte community homogenization between the 1940s and 2010s, likely due to modest land-use changes and minimal anthropogenic impacts. In Greece, previous studies have shown that hydromorphological and geomorphological features, such as bed substrate, habitat type, flow conditions, and channel slope, play a dominant role in structuring the taxonomic and functional composition of macrophyte communities, often outweighing the influence of physicochemical variables (Stefanidis et al., 2021a; 2023). Similarly, hydromorphological alterations have been reported to exert a stronger effect on the IBMR_GR than changes in water quality, highlighting the primary role of physical habitat conditions in shaping macrophyte-based ecological assessments (Stefanidis et al., 2022). We consider it highly likely that cumulative hydromorphological changes in lotic ecosystems of Greece (Stefanidis et al., 2020), often associated with channel simplification and habitat loss, could explain the observed temporal patterns in community composition. Moreover, recent climate change-driven increases in the frequency and intensity of extreme precipitation and flood events (Arnell and Gosling, 2016; Alifu et al., 2022), such as the flooding caused by Storm Daniel in 2023, which delivered record rainfall and exceptional river discharges across central Greece (Dimitriou et al., 2024; Kakavas et al., 2025), underscore how altered hydrological regimes can further disrupt river ecosystems and habitat structure, potentially influencing aquatic communities.

In summary, our results showed that earlier years (2014 and 2015) were compositionally distinct from later years (2021, 2022, and 2023), indicating substantial temporal turnover and potential directional changes in community composition. In contrast, consecutive recent years displayed lower and often non-significant dissimilarity, suggesting a relative stabilization of community structure during the most recent sampling period. However, studies investigating temporal patterns of β-diversity dynamics are particularly important but require large spatial and temporal scales (Magurran et al., 2019). In our case, although temporal variation in β-diversity was observed and a shift in 2022 was apparent, there is no strong evidence supporting a trend toward biotic homogenization. In contrast, our findings (e.g. LCBD variation) suggest a shift to more diverse communities. Yet, we have to point out that such analyses generally require long-term datasets to accurately detect directional changes, and a dataset spanning only 10 years may be insufficient to capture meaningful temporal trends. The sampling years included in this study span substantial intervals, and the sites sampled were not entirely consistent across all years. Although these gaps and inconsistencies reflect the practical realities of monitoring, they introduce limitations when interpreting long-term trends. Therefore, observed temporal patterns should be considered cautiously, as they may partly reflect the uneven sampling design rather than continuous ecological changes. To this end, we expect that the next monitoring program, which began in 2025, will provide new data and additions to the aquatic flora of Greek rivers, potentially revealing hidden patterns and trends.

Macrophyte monitoring within the context of the WFD presents significant challenges in terms of practical implementation, data collection and analysis, and result interpretation. Sampling biased toward indicator species may lead to underestimation of species richness and simplification of communities. Moreover, sites are often selected based on the intensity of anthropogenic stress that may be responsible for ecological degradation, making monitoring necessary. Such biases in plant and site selection can result in unreliable long-term data. However, the use of harmonized field methods targeting the entire aquatic flora can address this issue and ensure the collection of representative community samples.

Monitoring data can then be used to explore multiple aspects of biodiversity in addition to ecological quality assessments. Complementing traditional field studies, new tools such as eDNA monitoring could enhance the effectiveness of macrophyte surveys and improve our understanding of macrophyte composition at the scale of entire catchment (Suren et al., 2024; Espinosa Prieto et al., 2025). Ensuring consistent collection of macrophyte data is essential for long-term assessment of ecological quality and community dynamics, both for evaluating the effectiveness of applied mitigation measures and for identifying biotic responses over time. On top of that, implementing beta diversity metrics to biodiversity assessments can provide us with useful insights and vital information on the community responses to environmental change (Socolar et al., 2016; Heino et al., 2024). Identifying spatial and temporal patterns of compositional change can be particularly useful for quantifying ecological uniqueness and biotic homogenization processes that may occur because of environmental perturbation (Dubois et al., 2020; Luukkonen et al., 2024).

In our study, the temporal beta diversity analysis showed that low turnover between consecutive years may suggest that local assemblages are compositionally stable from year to year, suggesting resilience to short-term environmental variability. However, over longer temporal intervals we found a significant increase in dissimilarity that likely reflects the cumulative effects of environmental pressures such as hydrological alterations, habitat destruction, organic pollution and climatic trends that progressively modify habitat conditions and drive species replacement or loss.

From a management perspective, this finding underscores the importance of maintaining long-term monitoring programs to capture gradual ecological changes that may be overlooked by short-term assessments. Detecting broader community shifts requires data spanning several years and therefore, conservation and restoration strategies should focus on preserving habitat diversity and addressing multiple stressors that may cumulatively contribute to long-term compositional change. In this context, beta diversity provides a valuable metric for evaluating the temporal and spatial dynamics of aquatic vegetation and for guiding adaptive management of freshwater ecosystems.