Multiple-Facet Diversity Patterns of Aquatic Vegetation in Lakes along a Trophic Gradient

The EU Water Framework Directive foresees the ecological assessment of surface waters against identified pressures. Nutrient loading is the main pressure impairing the ecological quality of lake ecosystems, and aquatic macrophytes are considered good indicators of ecological response. In this study, we statistically assessed different aspects of aquatic plant (macrophyte) diversity in response to different trophic levels in Mediterranean lakes. We used 5690 relevés of aquatic vegetation, distributed over 305 transects, sampled in 18 freshwater lake ecosystems during 2013–2016. Our results show a significant decrease in taxonomic alpha diversity in lakes with a total phosphorus content above 100 μg/L. Syntaxonomic diversity followed the species richness pattern as well. Functional richness decreased along the trophic gradient, while functional dispersion was higher in lakes with high trophic levels. Taxonomic and functional beta partitioning presented changes in assembly processes leading to greater community homogeneity in lakes with higher trophic levels. In summary, we found no redundancy between taxonomic and functional diversity indices. These results provide novel insights into aquatic plant assembly processes of impacted freshwater lakes needed to forward conservation and restoration practices.


Introduction
The importance of freshwater ecosystems as biodiversity hotspots around the globe is well established [1,2]. However, human-induced changes have been impacting the integrity and biodiversity of these ecosystems through increasing nutrient, sediment and contaminant loads, water abstraction, and biotic exchanges [1,3,4]. Effective conservation and restoration measures to address biodiversity decline are needed, and, in this respect, increased efforts have been made to identify ecosystem diversity changes in order to allow for ecological predictions [5,6]. In Europe, the Water Framework Directive (WFD; [7]) necessitates the ecological assessment of the structure and functioning of surface water ecosystems, and therefore holistic ecosystem assessment approaches are required [8,9]. However, most WFD compliant monitoring and assessment methods focus on the taxonomic aspect of freshwater ecosystems [10]. Metrics taking into consideration all attributes of diversity, including its functional component, failed to find a place in most biomonitoring protocols of freshwater ecosystems [11][12][13]. It has been emphasized that in order to evaluate and quantify ecosystem functioning, new assessment tools should be developed [10].
Following Noss [14], who proposed the exploration of biodiversity's three main facets (compositional, structural, and functional), research started to focus on multiple-facet diversity approaches [15,16]. These approaches aim to evaluate biodiversity changes by using various diversity indices. Taxonomic diversity (TD), the local or regional species pool [15,17], constitutes the most commonly used. TD is further divided into taxonomic alpha diversity, i.e., species richness in a specific site, and taxonomic beta diversity, i.e., variation in species composition among sites [18,19]. Changes in taxonomic beta diversity reflect two different phenomena: one is the species turnover from one site to another, and the other is species absence when comparing sites, known as nestedness [17,18,20,21]. Another important component of biodiversity is functional diversity (FD), which suggests that the species are not equivalent regarding ecosystem processes [16,22]. Likewise, FD consists of an alpha component, which contains the functional richness in a site, and a beta one, which stands for the variation in functional trait composition among sites of a region, further divided into functional turnover (i.e., functional replacement) and functional nestedness (i.e., functional loss) [20,23,24]. Finally, another level of biodiversity is community or syntaxonomic diversity (SD), using plant communities (syntaxa) [25,26].
Aquatic macrophyte diversity, especially species richness, and its relation with environmental conditions in freshwater ecosystems has been previously examined in several lake ecosystems [27][28][29]). Studies proved that they respond well to eutrophication pressure [30][31][32], which had, as a consequence, to classify them as bioindicators in WFDcompliant assessment systems. However, there was not enough focus on multi-scale patterns and community diversity in relation to eutrophication pressure [11][12][13]33]. Furthermore, aquatic macrophyte diversity has not received proper attention in the Mediterranean region, and research gaps reflecting different monitoring traditions in southern, central and northern European countries, in the context of WFD, are still profound [11,[33][34][35][36].
In this study, we aimed to identify and describe different aspects of aquatic plant diversity in Mediterranean lakes in response to different trophic levels. More specifically, our first objective was to explore the patterns of the taxonomic, syntaxonomic and functional diversity of aquatic vegetation along the trophic gradient, by means of calculating a number of diversity indices for lakes in different trophic levels. Then, the resulting diversity patterns were used for our second objective, which was to infer the changes in aquatic plant assembly rules among the different lake trophic levels.

Study Area and Vegetation Data
Aquatic and littoral vegetation data were collected from 18 lakes ( Figure 1 and Table 1) belonging to the Greek National Water Monitoring Network (GNWMN) [37].  Table 1 for lake names.   Table 1 for lake names. All 18 lakes belong to Mediterranean lake types [39] and are distributed throughout Greece. They are subject to a wide range of pressures, including different nutrient loads and land use in the catchment areas [40,41]. Each lake was surveyed once in May to September 2013-2016 during the main growing season of aquatic and riparian plants (Table 1), by applying the belt transect-mapping method (for more details see [37]). We collected 5690 vegetation relevés of 4 m 2 , distributed among 305 transects along lake depth gradients (Table 1). All macrophytes (angiosperms, pteridophytes, bryophytes, charophytes and filamentous green algae) were recorded and determined to (sub-)species level (except filamentous green algae). Their abundance at plot level was estimated using the semi-quantitative five-point DAFOR scale [42], transformed to percentage cover values as follows: Dominant = 87.5%, Abundant = 50%, Frequent = 17.5%, Occasional = 5.5% and Rare = 0.5%.

Lake Trophic Status Parameters
Environmental data (e.g., seasonal total phosphorus concentration in the water column) were collected from each lake in the context of GNWMN (see Zervas et al. [41] for more details). Based on these data, we calculated annual mean total phosphorus (TP) concentration in the water column of each lake. We considered TP concentration as proxy to eutrophication pressure [43,44], and we used it in order to assess the relationships between aquatic vegetation and lake trophic level.
Various trophic status classification schemes have been developed in order to connect the eutrophication levels of water bodies and nutrient concentrations [45]. In the context of this study, the lakes were grouped into four categories according to their TP content as follows: low trophic group, with TP < 20 µg/L, which is considered a threshold value for reference lakes [43,46]; moderate trophic group, with TP 20-50 µg/L, which is the range used to derive the good/moderate boundaries for Mediterranean lakes in WFD assessment systems [47,48]; very high trophic group, TP > 100 µg/L, which responds to hyper-eutrophic lakes according to the OECD classification scheme [49]; and the in-between range, of the high trophic group, with TP 50-100 µg/L.

Diversity Analysis
Several diversity indices (i.e., taxonomic alpha and beta diversity, syntaxonomic diversity, functional alpha and beta diversity) were calculated for the aquatic and littoral vegetation recorded in each lake. For the purpose of visualizing the diversity patterns and turning points, these indices were plotted against the lake trophic gradient. In order to acquire diversity values with respect to different lake sizes, calculations of diversity indices (besides rarefaction curves) were performed on a transect level. This approach was also favored in order to overcome species number restrictions when calculating functional diversity indices [50,51]. Same-transect relevés were integrated, transect sizes were standardized in order to equalize the number of their relevés, and macrophyte taxa cover values were calculated for each lake transect.
For taxonomic alpha diversity, rarefaction curves and diversity profiles based on the accumulation of relevés (sampling-unit-based incidence data) for each lake were created with interpolation/extrapolation analysis in the iNEXT package for R [52] using Hill numbers, as proposed by Chao et al. [53]. Species richness, the Shannon diversity index, and the Simpson diversity index were also calculated for each transect using the vegan package for R [54].
In order to calculate syntaxonomic diversity, i.e., the richness of plant communities in each transect, we used the results of the phytosociological analysis in Zervas et al. [55]. Each relevé was coded according to its community affiliation, and syntaxa richness in each transect was calculated using vegan package for R [54].
For functional alpha diversity, we calculated the multidimensional Functional Richness index (FRic: amount of niche space occupied by taxa) of Villéger et al. [51] and Rao's quadratic entropy (FD Q : dispersion of taxa in trait niche space) [56], as proposed by Lal-iberté and Legendre [57]. Calculations were based on a multiple-trait dataset with traits covering life forms, growth morphology, survival characteristics, ecological preferences, reproduction and dispersal attributes (for more details, see Zervas et al. [41]). The trait dataset was transformed into a species distance matrix using Gower dissimilarity [58] for categorical and ordinal trait values, and was analyzed by principal coordinates analysis (PCoA) in order to reduce the trait dimensions to two axes. The resulting PCoA axes were used in calculating FRic and FD Q diversity indices for each transect of the study. All above-mentioned analysis was executed using the FD package for R [59].
Taxonomic beta diversity in each lake was computed by applying the multiple-site extension to pairwise partitioned beta diversity, as proposed by Ensing and Pither [60]. Abundance community data of all transects in each lake were analyzed by using the betapart package for R [61]. This process calculates the variation in species composition among those transects and partitions the variation in species replacement (turnover) and species loss (nestedness) from transect to transect in each lake.
Functional beta diversity in each lake was computed accordingly, by applying the multiple-site beta diversity extension on the abundance community data of all transects in each lake. For the required trait input, we used the same minimized dimension matrix we received from the PCoA analysis of the trait data in the functional diversity calculations. Transect functional beta diversity values, and their partitioning to functional turnover and functional nestedness, were calculated by using the betapart package for R [61].
For each diversity index, one-way ANOVA (analysis of variance, [62]) was applied on the average diversity values of transects belonging to lakes in the different trophic groups, in order to evaluate the diversity differences along the trophic gradient. The Tukey Honest Significant Differences post-hoc test (TukeyHSD; [63]), with a 95% family-wise confidence level, was used to explore significant differences of diversity values between different lake trophic levels. ANOVAs and TukeyHSD tests were made using the stats core-package for R [64]. Furthermore, correlations among all diversity metrics and TP were assessed in order to (i) investigate the strength of the relation of each diversity metric independently with TP as an indicator of lake trophic level and (ii) explore the redundancy in calculating different diversity metrics. Taking into consideration the nonlinearity of diversity responses to environmental gradients, a non-linear correlation analysis [65] was applied among all the estimated transect diversity values and the lake TP values by using the nlcor function for R [64].
All calculations were performed in R environment version 3.5.3 [64].

Results
Regarding the lake classification in trophic groups, six out of the 18 lakes (Amvrakia, Kourna, Trichonida, Feneos, Paralimni, and Yliki) were grouped into the low trophic group (<20 µg/L TP). Another six lakes (Doirani, Megali Prespa, Ozeros, Kastoria, Mikri Prespa, and Vegoritida) were assigned to the moderate trophic level (20-50 µg/L TP). Three lakes (Chimaditida, Petres, and Volvi) were placed in the high trophic group (50-100 µg/L TP), and the remaining three (Lysimachia, Pamvotida, and Zazari) in the very high trophic group (>100 µg/L TP). The mean depth of low trophic group lakes fluctuated between 5 and 29 m, between 4 and 26 m for moderate trophic group lakes, between 1 and 13 m for high trophic group lakes, and between 4 and 5 m for very high trophic group lakes.
Rarefaction curves ( Figure 2) provided a diversity profile for each lake, as well as the four trophic groups they belong to. Diversity indices for most lakes reached a plateau during the interpolation of their relevés sampled, except for Yliki (low trophic) in all three indices, and Paralimni (low trophic), Lysimachia (very high trophic), Pamvotida (very high trophic) and Zazari (very high trophic) for species richness. Variance in the diversity values of lakes in the low and moderate trophic groups was greater than that observed in lakes in the high and very high trophic groups. Paralimni in the low trophic group was the most species-rich lake, followed by Megali Prespa and Doirani in the moderate trophic group, and Trichonida also in the low trophic one. Despite the great variation, lakes belonging to low and moderate trophic groups combined had higher species richness than those in the high and very high trophic groups. Shannon and Simpson diversity indices showed a clearer image of the diversity differences between the lakes representing different trophic groups. Doirani, Vegoritida and Megali Prespa of the moderate trophic group were the most diverse lakes, followed by lakes in the low and high trophic groups. Lakes belonging in the very high trophic group presented the lowest diversity values.
indices, and Paralimni (low trophic), Lysimachia (very high trophic), Pamvotida (very high trophic) and Zazari (very high trophic) for species richness. Variance in the diversity values of lakes in the low and moderate trophic groups was greater than that observed in lakes in the high and very high trophic groups. Paralimni in the low trophic group was the most species-rich lake, followed by Megali Prespa and Doirani in the moderate trophic group, and Trichonida also in the low trophic one. Despite the great variation, lakes belonging to low and moderate trophic groups combined had higher species richness than those in the high and very high trophic groups. Shannon and Simpson diversity indices showed a clearer image of the diversity differences between the lakes representing different trophic groups. Doirani, Vegoritida and Megali Prespa of the moderate trophic group were the most diverse lakes, followed by lakes in the low and high trophic groups. Lakes belonging in the very high trophic group presented the lowest diversity values.  Table 1 for lake codes) divided into four trophic groups (low: low trophic group; mod: moderate trophic group; high: high trophic group; v. high: very high trophic group). The solid part of the lines predicts the rate of accumulation of the diversity metric values in each lake with increasing number of sampling units (relevés) by interpolation of data, while the dashed part of the lines continues the same predictions by extrapolation. The curve of the Yliki lake (low trophic group) is hardly visible because of its small length (small number of relevés) and its overlap with the curves of other lakes.
Alpha diversity patterns along the lake trophic gradient are summarized in Figure 3. Each bar represents the mean diversity value of all transects belonging to lakes within the same trophic group. One-way ANOVA showed that there were statistically significant (p < 0.001) differences among all three alpha diversity values of transects belonging to the different trophic lake groups. Greater values of species richness were calculated in lakes with a moderate trophic level, while greater values of Shannon and Simpson diversity values were calculated in lakes with a high trophic level. Lower values for all three alpha diversity metrics were found in very high trophic level lakes. However, the TukeyHSD test showed significant differences only between the lakes in the very high trophic group  Table 1 for lake codes) divided into four trophic groups (low: low trophic group; mod: moderate trophic group; high: high trophic group; v. high: very high trophic group). The solid part of the lines predicts the rate of accumulation of the diversity metric values in each lake with increasing number of sampling units (relevés) by interpolation of data, while the dashed part of the lines continues the same predictions by extrapolation. The curve of the Yliki lake (low trophic group) is hardly visible because of its small length (small number of relevés) and its overlap with the curves of other lakes.
Alpha diversity patterns along the lake trophic gradient are summarized in Figure 3. Each bar represents the mean diversity value of all transects belonging to lakes within the same trophic group. One-way ANOVA showed that there were statistically significant (p < 0.001) differences among all three alpha diversity values of transects belonging to the different trophic lake groups. Greater values of species richness were calculated in lakes with a moderate trophic level, while greater values of Shannon and Simpson diversity values were calculated in lakes with a high trophic level. Lower values for all three alpha diversity metrics were found in very high trophic level lakes. However, the TukeyHSD test showed significant differences only between the lakes in the very high trophic group and the rest of them for the species richness and Shannon indices, while it showed significant differences between the lakes in the very high trophic group and those in the high trophic group, as well as those in the moderate and low trophic groups combined for the Simpson index. and the rest of them for the species richness and Shannon indices, while it showed significant differences between the lakes in the very high trophic group and those in the high trophic group, as well as those in the moderate and low trophic groups combined for the Simpson index. and Rao's quadratic entropy in lake transects, grouped according to annual mean total phosphorus (TP) content in each lake's water column (low, moderate, high and very high trophic groups). Standard error bars for all mean values are given. One-way ANOVA results (degrees of freedom, F statistic, and p-value) for each diversity metric are presented in each graph's headline. Statistically significant differences of diversity means in consecutive trophic levels, according to Tuk-eyHSD multiple comparison for 95% family-wise confidence level, are indicated by different alphabet letters among a, b, c and d.
Syntaxonomic diversity (number of distinct plant communities) followed the same pattern with species richness (Figure 3), peaking in the moderate trophic group. One-way ANOVA showed that there was also a statistically significant (p < 0.001) difference among the values belonging to lakes of different trophic levels, but TukeyHSD differentiated statistically only the very high trophic group from the others.
Functional diversity indices, Functional Richness and Rao's quadratic entropy revealed a different pattern (Figure 3). According to one-way ANOVA, there were statistically significant (p < 0.001) differences among the functional diversity values of lakes belonging to different trophic groups. Significant higher values of functional richness were calculated for lakes belonging to the low trophic group, moderate values were calculated for moderate and high trophic groups, and lower values were calculated for those in the very high trophic group. The TukeyHSD test verified this differentiation among the trophic groups (low, moderate and high combined, and very high groups). On the other and Rao's quadratic entropy in lake transects, grouped according to annual mean total phosphorus (TP) content in each lake's water column (low, moderate, high and very high trophic groups). Standard error bars for all mean values are given. One-way ANOVA results (degrees of freedom, F statistic, and p-value) for each diversity metric are presented in each graph's headline. Statistically significant differences of diversity means in consecutive trophic levels, according to TukeyHSD multiple comparison for 95% family-wise confidence level, are indicated by different alphabet letters among a, b, c and d.
Syntaxonomic diversity (number of distinct plant communities) followed the same pattern with species richness (Figure 3), peaking in the moderate trophic group. One-way ANOVA showed that there was also a statistically significant (p < 0.001) difference among the values belonging to lakes of different trophic levels, but TukeyHSD differentiated statistically only the very high trophic group from the others.
Functional diversity indices, Functional Richness and Rao's quadratic entropy revealed a different pattern (Figure 3). According to one-way ANOVA, there were statistically significant (p < 0.001) differences among the functional diversity values of lakes belonging to different trophic groups. Significant higher values of functional richness were calculated for lakes belonging to the low trophic group, moderate values were calculated for moderate and high trophic groups, and lower values were calculated for those in the very high trophic group. The TukeyHSD test verified this differentiation among the trophic groups (low, moderate and high combined, and very high groups). On the other hand, Rao's quadratic entropy was found to receive higher values for lakes in the high trophic group and lower ones for the lakes in the low and very high trophic groups. The TukeyHSD test revealed significant differences in Rao's values for all trophic groups. The summary of the results regarding taxonomic and functional beta diversity and their partitions, turnover and nestedness, are shown in Figure 4. Both one-way ANOVAs for taxonomic and functional beta values showed that there were statistically significant (p < 0.001) differences among the mean beta diversity values of lakes belonging to different trophic groups. For taxonomic beta diversity, TukeyHSD showed that differences between low and moderate trophic groups, and moderate and high-very high trophic groups were statistically significant. For functional beta diversity, TukeyHSD showed significant differences among all consecutive trophic groups. Beta partitioning showed more clear patterns along the trophic gradient. Taxonomic beta turnover received higher values in lakes with moderate and high trophic levels, and lower ones in those of low and very high trophic levels, while taxonomic beta nestedness followed the reverse pattern. For both partitions, ANOVA showed statistically significant (p < 0.001) differences among the mean diversity values, and TukeyHSD recognized significant differences between the values of low and moderate-high trophic groups and between those of moderate-high and very high trophic groups. Moreover, taxonomic beta turnover values were higher (>0.5) than those of taxonomic beta nestedness (<0.3) for all trophic groups. Functional beta turnover was significantly lower in the high trophic group, where functional beta nestedness peaked. ANOVA showed statistically significant (p < 0.001) differences among the mean index values for both functional beta partitions along the trophic gradient, and TukeyHSD recognized significant differences among the diversity values in lakes with low-moderate trophic levels together, as compared to those in the high trophic group, and to those in the very high trophic group. Moreover, functional beta turnover values were found to be lower (<0.35) than those of functional beta nestedness (>0.35) for all trophic groups. hand, Rao's quadratic entropy was found to receive higher values for lakes in the high trophic group and lower ones for the lakes in the low and very high trophic groups. The TukeyHSD test revealed significant differences in Rao's values for all trophic groups. The summary of the results regarding taxonomic and functional beta diversity and their partitions, turnover and nestedness, are shown in Figure 4. Both one-way ANOVAs for taxonomic and functional beta values showed that there were statistically significant (p < 0.001) differences among the mean beta diversity values of lakes belonging to different trophic groups. For taxonomic beta diversity, TukeyHSD showed that differences between low and moderate trophic groups, and moderate and high-very high trophic groups were statistically significant. For functional beta diversity, TukeyHSD showed significant differences among all consecutive trophic groups. Beta partitioning showed more clear patterns along the trophic gradient. Taxonomic beta turnover received higher values in lakes with moderate and high trophic levels, and lower ones in those of low and very high trophic levels, while taxonomic beta nestedness followed the reverse pattern. For both partitions, ANOVA showed statistically significant (p < 0.001) differences among the mean diversity values, and TukeyHSD recognized significant differences between the values of low and moderate-high trophic groups and between those of moderate-high and very high trophic groups. Moreover, taxonomic beta turnover values were higher (>0.5) than those of taxonomic beta nestedness (<0.3) for all trophic groups. Functional beta turnover was significantly lower in the high trophic group, where functional beta nestedness peaked. ANOVA showed statistically significant (p < 0.001) differences among the mean index values for both functional beta partitions along the trophic gradient, and TukeyHSD recognized significant differences among the diversity values in lakes with low-moderate trophic levels together, as compared to those in the high trophic group, and to those in the very high trophic group. Moreover, functional beta turnover values were found to be lower (<0.35) than those of functional beta nestedness (>0.35) for all trophic groups. With respect to the relationships between all diversity metrics and TP (Table 2), all paired-metric correlations were statistically significant according to non-linear correlation analysis (adjusted p < 0.05). However, low (<0.8) correlation coefficients were calculated for most relationships. None of the individual diversity metrics showed a high correlation with TP. Between the different diversity metrics, species richness was highly correlated with the Shannon index and syntaxonomic richness, while the Shannon index was also highly correlated with the Simpson index and syntaxonomic richness. Furthermore, high correlations were also found between the different beta diversity metrics. On the other hand, no high correlations were found between taxonomic and functional diversity metrics. Table 2. Overview of the non-linear correlation analysis results among all diversity metrics (SpeRich: Species richness; Shan: Shannon index; Simp: Simpson index; Syntax: Syntaxonomic richness; FunRich: Functional richness; Rao: Rao's quadratic entropy; BetOver: Overall taxonomic beta diversity index; BetTurn: Taxonomic beta turnover; BetNest: Taxonomic beta nestedness; FunBeta: Overall functional beta diversity index; FunTurn: Functional beta turnover; FunNest: Functional beta nestedness) and total phosphorus content (TP) in lakes. All paired-metric correlations were found statistically significant (adjusted p < 0.05).

Taxonomic Alpha Diversity Patterns
Grime [66] described the "hump-back" model of diversity, in which high biodiversity is observed at intermediate intensities of disturbance due to the retardation of competitive exclusion. This "hump-back" model was also discovered governing the relationship between aquatic plant species richness and trophic level, recording higher species richness in mesotrophic lakes, and lower in nutrient-poor and nutrient-rich ones [27,29,67]. In this intermediate range of nutrient loads, an increase in emergent macrophyte richness and abundance due to lower levels of competitive exclusion overcompensates for the decrease in submerged macrophyte richness and abundance due to elevated water turbidity [67][68][69]. Jeppesen et al. [70] placed this intermediate range of nutrient loads between 50 and 150 µg/L TP for northern temperate lakes, while Romo et al. [71] suggested 50 µg/L TP as a critical threshold for Mediterranean lakes to avoid the shift to a turbid-water phase and preserve macrophyte dominance.
Our results indicate that there is a similar but less obvious pattern for the taxonomic alpha diversity of aquatic plant species in Mediterranean lakes, with higher species richness in lakes with moderate trophic levels (20-50 µg/L TP) and higher Shannon and Simpson index values for lakes with high trophic levels (50-100 µg/L TP). However, only the Simpson index was able to reveal a statistically significant peak of diversity in lakes with a high trophic level, while species richness and the Shannon index revealed significant diversity changes only between lakes below and above the 100 µg/L TP threshold (highvery high lake trophic boundary). The "hump-back" pattern might have been more profound in our study if data concerning very low trophic level lakes (<10 µg/L TP), which are commonly used in other studies (e.g., [30,69]), were available in our dataset. Parts of the "hump-back" diversity pattern were described recently by García-Girón et al. [36] for macrophytes in Mediterranean ponds, as well as in studies of other aquatic organisms, such as phytoplankton, zooplankton, fish and macroinvertebrates [72][73][74].
Penning et al. [30] and Kolada et al. [69] were also able to find the same "hump-back" response, with macrophyte species richness peaking at 20-50 µg/L TP, when testing its indicator value for a great number of northern and central European lakes. However, great variation in their dataset and the resulting weak relation to TP content led to the conclusion that species richness alone is of limited explanatory power. Similarly, we found great variation in the values of species richness, Shannon and Simpson indices, especially for low and moderate trophic level lakes, leading to uncertainty in the assessment results. Moreover, neither of these three alpha diversity metrics independently showed a high correlation with TP lake content. Similar results presented by García-Girón et al. [36] strengthened their suggestions for the complementary use of several diversity metrics. Diversity metrics for aquatic vegetation are considered sensitive to environmental "noise" associated with a number of different environmental drivers related to habitat conditions, habitat heterogeneity and other water physico-chemical characteristics [29,[33][34][35]. Therefore, our results concerning alpha diversity metrics can be attributed to two different scenarios. The first scenario is that perhaps nutrient loading was not the most important driver of the spatial variation of macrophyte communities within a lake, and other drivers such as habitat heterogeneity or other local disturbance factors should be taken into consideration [75,76]. The second scenario is that nutrient loading could be the most important driver of macrophyte spatial variation but it does not affect the whole lake body homogenously, having in mind that many lakes may be impacted by point-source nutrient loadings and diverse water quality conditions prevail among distant shoreline habitats [76][77][78].

Taxonomic Beta Diversity Patterns
Another important aspect when analyzing diversity patterns for conservation purposes is the spatial variation of species composition among different communities, known as taxonomic beta diversity [18,19,79]. Its two components, species turnover (or species replacement) and nestedness (species loss), can reflect the drivers behind spatial variation, such as the natural and life history of species (e.g., distribution dynamics, dispersal strategies and evolution), species interactions (e.g., competition, mutualism and parasitism), and other ecological factors (e.g., environmental conditions, habitat heterogeneity and disturbance) [79][80][81].
According to our results, the spatial variation of species composition (total beta diversity), calculated among multiple sites within each lake, remains high at all trophic levels, following the results of similar studies [36,82,83]). Partitioning beta diversity revealed two distinct patterns of species turnover and nestedness along the trophic gradient. Species turnover among sites in lakes with moderate and high trophic levels was higher than that in lakes with low and very high trophic levels, where nestedness increased. This indicates that spatial variation in moderate and high trophic level lakes is attributed mainly to species replacement from site to site, maintaining high species richness in those lakes. On the other hand, spatial variation in lakes with a low trophic level with slightly lower species richness, and in lakes with a very high trophic level with significantly lower species richness, is attributed to a mix of species replacement and species loss from site to site. Taking into consideration that species turnover is often connected with environmental filtering and competition, while nestedness indicates other ecological processes, such as human disturbance [21,79,81,84], we conclude that ecological factors other than nutrient enrichment are also playing a significant role in macrophyte spatial distribution patterns within lakes (scenario 1), and local site nutrient profiles are important (scenario 2) in order to disentangle and quantify the relative contribution of each ecological process to these spatial distribution patterns.

Syntaxonomic Diversity Patterns
Syntaxonomic (or community) diversity is another way to express vegetation diversity by presenting at the same time plant diversity and their specific habitat features, thus acting as an indicator for ecological heterogeneity [25,26]. Plant species are the basis for describing syntaxonomic (phytosociological) units; therefore, syntaxonomic richness may be greatly affected by species richness. However, the number of plant communities in a habitat does not always imply the same number of species (species richness), since plant species are not equally distributed in communities, and disturbance may result in an increase/decrease in companion taxa without affecting the number of their communities [25,26].
In our results, however, we found that syntaxonomic richness followed a similar pattern to that of species richness (highly correlated), i.e., slightly higher numbers of plant communities were observed in lakes with moderate trophic levels, slightly lower numbers in lakes with low and high trophic levels, and significantly lower numbers in lakes with very high trophic levels. This result may be attributed to the fact that macrophyte communities contain a small number of taxa [55] and reinforces the "hump-back" model of diversity of macrophyte communities in lakes in the middle of the trophic gradient. However, variation in syntaxonomic richness values was also an issue for lakes with low and moderate trophic levels, reinforcing the need for the identification of parameters other than nutrient loading.

Functional Alpha Diversity Patterns
It has been extensively discussed that ecosystem functioning often depends on the richness and distribution of species, taking into consideration their functional role in the ecosystem [6,85,86]). High species richness does not always imply high ecosystem resilience (e.g., nuisance species) and it is not always the most desired target for conservation and restoration purposes [5,87,88]. Functional diversity is another facet of biodiversity adopted in order to evaluate diversity responses to ecosystem disturbance [87][88][89].
A number of different indices have emerged in order to quantify functional diversity, by assessing the value and range of organismal traits in ecosystems [50,51,90]. Functional richness, i.e., the functional niche space occupied by a species assemblage, and functional divergence, i.e., the distance of high species abundances from the center of the functional niche space, have been suggested as the best solution in analyzing community assembly processes under stress gradients [90]. Functional richness is increasing when present species occupy a larger functional niche space [50,51,90]. Therefore, functional richness is usually positively correlated with species richness (having more species leads to a bigger functional space) but not always, as the functional traits of some community assemblies may be more closely clustered than others [50]. Functional divergence is increasing when the majority of the community (species present and their abundance) is occupying the edges of the functional niche space [50,90,91]. Rao's quadratic entropy is an index that measures both functional richness and divergence (conceptually similar to functional dispersion), increasing when either the functional niche space occupied by species increases and/or when these species are displaced further from the functional niche center [57,90,91].
In our results, we witnessed a drop in functional richness in the lakes of the moderate and high trophic groups, in relation to those in the low trophic group, while maintaining similar levels of species richness. This indicates that macrophyte species found in low trophic level lakes are being replaced in moderate and high trophic level ones by species that occupy a smaller functional niche space; thus, such species are defined by more homogenous trait patterns. Lakes with very high trophic levels present simultaneously lower functional and lower species richness, indicating both species and functional loss. This result of ours only partly follows the findings of previous studies [92][93][94]) that suggested species loss as the primal driver for functional changes along trophic gradients. On the other hand, Rao's quadratic entropy revealed the similar "hump-back" response, much like the taxonomic indices, concurring with the results of García-Girón et al. [36]. This indicates that functional richness in low trophic level lakes is inflated by rare species with extreme trait values, while the majority of the community occupies a much smaller functional niche space. The highest Rao index values in high trophic level lakes show that their macrophyte communities are more dispersed in the functional niche space, following the evenness pattern presented by the Simpson index. Therefore, while low trophic level lakes present a greater amount of functional richness, they are characterized by lower functional differentiation than needed to ensure ecosystem resilience to changing environmental conditions [88,89,95].

Functional Beta Diversity Patterns
Much like taxonomic beta diversity, the spatial variation of functional trait compositions can be analyzed by functional beta diversity and its two analogous components, functional turnover (i.e., functional replacement) and functional nestedness (i.e., functional loss) [20,23,24]. Functional beta diversity metrics attempt to detect patterns in niche-based assembly processes, reflecting the responses of traits rather than species to the same ecological drivers that taxonomic beta diversity attempts to investigate [20,96].
According to our results, there were changes in spatial variation in trait composition (total functional beta diversity) along the trophic gradient, but we were not able to detect the profound downgrading pattern presented by Fu et al. [97]. The same factors that were suggested for taxonomic beta diversity, i.e., factors other than nutrient loading environmental drivers of macrophyte assemblages' variation (scenario 1), and differentiated local trophic conditions inside the lakes (scenario 2), should be taken into consideration in order to quantify beta diversities. Partitioning functional beta diversity showed that there is an increase in spatial variation due to functional loss in high trophic level lakes, indicating that species-rich macrophyte communities in some sites (transects) of those lakes are losing parts of their functional niche space. Concurring with the results of Fu et al. [97], the overall functional beta diversity in all trophic levels was primarily attributed to the functional nestedness component (functional loss), while taxonomic beta diversity was mostly attributed to species turnover (species replacement). This indicates that along the trophic gradient, regardless of changes in species richness, species replacement among sites is the dominating factor of spatial variation, accompanied at the same time by functional loss. This concept agrees with a number of previous studies [98][99][100]) suggesting that nutrient loads lead to increased homogeneity in freshwater biotic assemblages.

Conclusions
The response of different aspects of aquatic plant diversity along a trophic gradient found in 18 Mediterranean lakes was assessed. A number of different taxonomic, syntaxonomic, and functional indices were applied and their patterns of change along different trophic levels were investigated. Higher values of species richness and the Shannon and Simpson diversity indices were recorded in lakes with TP levels up to 100 µg/L TP, and significantly lower values were recorded in lakes with TP levels higher than 100 µg/L. The "hump-back" diversity pattern was identified statistically only for the Simpson diversity index, perhaps due to the lack of lakes with very low trophic levels. Syntaxonomic richness followed the species richness pattern. Great variations in diversity values were found probably due to environmental drivers other than nutrient loading. Functional richness followed a decreasing pattern along the trophic gradient, while Rao's index showed a "hump-back" pattern peaking in lakes with high trophic levels, indicating that low trophic level lakes contain rare species with extreme trait values, while the majority of the community abundance occupies a much smaller functional niche space. The partitions of taxonomic and functional beta diversity, turnover and nestedness described changes in assembly processes leading to greater community homogeneity in lakes with higher trophic levels. None of the diversity indices independently were highly correlated with lake trophic content, and no high correlation (redundancy) was found among taxonomic alpha or beta diversity indices and functional alpha or beta diversity indices. In conclusion, the complementary use of several diversity metrics is recommended in order to identify changes in community assembly processes as a response to eutrophication pressure.

Data Availability Statement:
The data that support the reported results of this study were used under license from The Goulandris Natural History Museum, Greek Biotope/Wetland Centre (EKBY). They are available upon reasonable request from DZ and with permission of The Goulandris Natural History Museum, Greek Biotope/Wetland Centre (EKBY).