Responses of Freshwater Diatoms and Macrophytes Rely on the Stressor Gradient Length across the River Systems

: Phytobenthic diatoms and macrophyte communities respond differently to stressors in aquatic environments. For the assessment of the ecological status of rivers in Slovenia, we use several indices, including the River Macrophyte Index (RMI) and Trophic index (TI) based on macrophyte and phytobenthic diatoms communities, respectively. In the present study, we examined the relationships between nutrient variables and values of RMI and TI using varied stressor gradient lengths. We also aimed to explain the variability of macrophyte and diatom communities with different stressors, namely nutrients and land cover variables and their combinations. The relationships of RMI and TI with nutrient variables varied signiﬁcantly and were affected by the length of the stressor gradient. We obtained a stronger relationship between the RMI and total phosphorous at an approximately <0.3-mg/L annual mean value, while, for the relationships with the TI, the values were signiﬁcant at bigger gradient lengths. The greatest share of variability in the macrophyte and diatom community was explained by the combination of land use and nutrient variables and the lowest share by phosphorus and nitrogen variables. When we applied a composite stressor gradient, it explained a similar share of the variability of both macrophyte and diatom communities (up to 26%). A principal component analysis (PCA) based on land use and nutrient stressor gradient revealed that the relationship between RMI EQR and PCA1 that represents intensive agriculture depends on the length of the gradient. The relationship was stronger for shorter gradients at lower values and decreased as the gradient extended towards higher values. Both tested assessment methods showed that macrophyte communities are more sensitive to shorter stressor gradients of lower values, whereas diatom communities are more sensitive to longer stressor gradient and higher values of the stressor.


Introduction
Phytobenthos and macrophytes are primary producers and, thus, at the base of the food web in rivers. Phytobenthos is a crucial element of the periphytic community that colonises different submerged substrates in the photic zone of the majority of aquatic ecosystems, including rivers. The diversity of the phytobenthos community is directly related to the heterogeneity of the habitats in the water body [1,2]. Their occurrence and abundance also depend on other factors, such as the availability of light and nutrients, water temperature, changes in water levels, and water velocity [3,4]. Human activity can significantly alter these parameters, which can consequently affect the structure and function of phytobenthic communities [5]. After disturbance, the colonization of the substrate with diatoms, which are important elements of phytobenthos, occurs relatively quickly, while the complex community develops only over a longer period [6,7]. Phytobenthos has The empirical research revealed relatively strong correlations between the nutrients and trophic indices based on phytobenthos; however, the correlations between nutrients and trophic indices based on macrophytes are rather weak [20,44,46,47]. This is possibly related to the length of the life cycle and growth form of both groups of primary producers. For macrophytes with a long life cycle and complex growth forms, the instantaneous chemical status of water could not be directly translated to the water quality during the macrophyte vegetation season. In addition, macrophytes are a very variable group regarding their growth forms, which also affects their response to various environmental factors [49][50][51], including water chemistry [52,53].
The aims of this study were (i) to examine the relationships between the nutrient variables and values of the River Macrophyte Index (RMI) and Trophic index (TI) using varied stressor gradient lengths and (ii) to explain the variability of the macrophyte and diatom communities with different stressors, namely nutrients and land cover variables and their combinations. We hypothesized (i) that TI will be related to nutrient concentrations in the water, while this will not be the case for RMI, (ii) that the relationships between stressor gradients and biotic indices will be affected by the length of the considered stressor gradients, and (iii) that a composite stressor gradient will explain a similar share of the variability of the macrophyte and phytobenthic communities.

Study Area and Sampling
Slovenia covers a total area of 20,273 km 2 and features 4573 km of rivers with catchments larger than 10 km 2 . The rivers extend over four European ecoregions: Po lowland (3), Alps (4), Dinaric Western Balkan (5), and Pannonian lowland (11) [54,55]. However, macrophytes in rivers are common and abundant especially in the ecoregions Dinaric Western Balkan, Pannonian lowland, and Po lowland. Our experimental set comprised macrophytes and phytobenthic diatoms data of 70 samples, obtained at 42 sampling sites ( Figure 1). Twenty-five sample sites were sampled once, 12 sample sites two times, 2 sample sites three times, and 3 sample sites were sampled five times. Half of the data were from Pannonian lowland, 32 from Dinaric Western Balkan, and 3 from Po Lowland. The sampling sites covered (near) natural-to-highly distorted conditions reflecting the various levels of perturbance caused by nutrients (phosphorous and nitrogen) and catchment land use.
Water 2021, 13, x FOR PEER REVIEW 3 of 21 annual measurements of the nutrients, while the nutrient content in the water can vary greatly during the year [48]. The empirical research revealed relatively strong correlations between the nutrients and trophic indices based on phytobenthos; however, the correlations between nutrients and trophic indices based on macrophytes are rather weak [20,44,46,47]. This is possibly related to the length of the life cycle and growth form of both groups of primary producers. For macrophytes with a long life cycle and complex growth forms, the instantaneous chemical status of water could not be directly translated to the water quality during the macrophyte vegetation season. In addition, macrophytes are a very variable group regarding their growth forms, which also affects their response to various environmental factors [49][50][51], including water chemistry [52,53]. The aims of this study were (i) to examine the relationships between the nutrient variables and values of the River Macrophyte Index (RMI) and Trophic index (TI) using varied stressor gradient lengths and (ii) to explain the variability of the macrophyte and diatom communities with different stressors, namely nutrients and land cover variables and their combinations. We hypothesized (i) that TI will be related to nutrient concentrations in the water, while this will not be the case for RMI, (ii) that the relationships between stressor gradients and biotic indices will be affected by the length of the considered stressor gradients, and (iii) that a composite stressor gradient will explain a similar share of the variability of the macrophyte and phytobenthic communities.

Study Area and Sampling
Slovenia covers a total area of 20,273 km 2 and features 4573 km of rivers with catchments larger than 10 km 2 . The rivers extend over four European ecoregions: Po lowland (3), Alps (4), Dinaric Western Balkan (5), and Pannonian lowland (11) [54,55]. However, macrophytes in rivers are common and abundant especially in the ecoregions Dinaric Western Balkan, Pannonian lowland, and Po lowland. Our experimental set comprised macrophytes and phytobenthic diatoms data of 70 samples, obtained at 42 sampling sites ( Figure 1). Twenty-five sample sites were sampled once, 12 sample sites two times, 2 sample sites three times, and 3 sample sites were sampled five times. Half of the data were from Pannonian lowland, 32 from Dinaric Western Balkan, and 3 from Po Lowland. The sampling sites covered (near) natural-to-highly distorted conditions reflecting the various levels of perturbance caused by nutrients (phosphorous and nitrogen) and catchment land use.   The biological and nutrient (nitrogen and phosphorous) data originated from a national surface water monitoring programme in Slovenia conducted by the Slovenian Environment Agency between 2007 and 2016. Macrophytes (vascular plants, bryophytes, and charophytes) were sampled during the main vegetation period (from June to September). The macrophyte sampling procedure followed the standardised Slovenian river bioassessment protocol [56]. Sampling was carried out on a 100-m-long river stretch during the low-to-medium water levels. Depending on the sampling site characteristics, sampling was performed by wading across the channel, from the bank of the stream, or from a boat. Submerged, floating, and emergent plants were recorded. Species abundance was estimated as the relative plant biomass using a five-degree scale: 1 = very rare, 2 = rare, 3 = common, 4 = frequent, and 5 = abundant and predominant [57,58]. Macrophyte species presence and abundance were used to calculate the River Macrophyte Index (RMI) [22], which is used for the assessment of the river ecological status and classification of river water bodies in Slovenia [37,59]. The RMI was calculated using the following equation [22]: where: Q Ai = abundance of the taxon i from the group A; Q ABi = abundance of the taxon i from the group AB; Q BCi = abundance of the taxon i from the group BC; Q Ci = abundance of the taxon i from the group C; Q Si = abundance of the taxon i from the groups A, AB, B, BC, and C; n A = total number of taxa in group A; n AB = total number of taxa in group AB; n BC = total number of taxa in group BC; n C = total number of taxa in group C; n S = total number of taxa in groups A, AB, B, BC, and C.
The phytobenthos sampling procedure followed the standardised Slovenian river bioassessment protocol [56]. The sampling procedure is based on the micro-habitat sampling approach in which one sample consists of subsamples from different habitats defined as a combination of the substrate, flow type, water depth, and riparian vegetation shading. In accordance with this methodology, phytobenthos subsamples were scraped with a sharp subject from various living and non-living solid underwater surfaces (pebbles, stones, rocks, macrophytes, submerged wood, etc.) that were at least 1-m distant from the bank. Collected subsamples were transferred into a bottle and preserved in 4% formaldehyde or alcohol. The determination of phytobenthos was out carried in the laboratory, where permanent slides were prepared, and then, at least 500 phytobenthic diatoms per sample were counted and identified. The presence and abundance of diatoms was used as a proxy for the phytobenthos community. These data were then used to calculate the Trophic index (TI) [46], which was used to assess the trophic status of rivers in Slovenia [59,60]. The TI was calculated using the following equation [46]: Gi = indicative weight of the taxon i; Hi = abundance of the taxon i.
Water samples for the nitrogen and phosphorous variables analyses were obtained monthly or at least four times a year (each season). In our analyses, only those variables were considered where data were available for all the selected sites: total nitrogen (TN), ammonium (NH 4 ), nitrate (NO 3 ), nitrite (NO 2 ), total phosphorus (TP), and orthophosphate (PO 4 3− ). For each variable, three variables were calculated: annual mean, annual maximum, and summer mean (Table 1). Land use characteristics were obtained at the catchment and subcatchment scales. The subcatchment areas were defined as parts of the catchment area upstream of a sampling site to the confluence of an important tributary (for details, see Pavlin [61]). Eight land use variables were calculated on the catchment and subcatchment scales using Corine Land Cover (CLC) data [62]. We distinguished four categories of land use: natural and semi-natural areas (CLC classes 3, 4, and 5); urban areas (CLC class 1); areas characterised by extensive agricultural practices (

Data Analyses
At the first step, we determined the relationship between the nutrient variables and biotic indices: the ecological quality ratio (EQR) values of the River Macrophyte Index (RMI EQR) and the Trophic Index (TI EQR) using the whole dataset (n = 70) (for details on the analytical procedure, see Figure 2). All relationships were defined using a Spearman rank correlation coefficient. Further on, we examined the influence of the stressor gradient length expressed as the mean annual total phosphorous range ( Figure 3) on the relationships between selected nutrient variables (ammonium maximum annual value and orthophosphate maximum annual value) and the biotic indices (RMI EQR and TI EQR). We determined the relationships using six different gradient lengths (Table 2). Due to the asymmetric distribution of the data, we determined the relationships based on the Spearman rank correlation. We used SPSS Statistics version 21.0(Armonk, NY, USA) to perform the analyses [63].
on the analytical procedure, see Figure 2). All relationships were defined using a Spearman rank correlation coefficient. Further on, we examined the influence of the stressor gradient length expressed as the mean annual total phosphorous range ( Figure 3) on the relationships between selected nutrient variables (ammonium maximum annual value and orthophosphate maximum annual value) and the biotic indices (RMI EQR and TI EQR). We determined the relationships using six different gradient lengths (Table 2). Due to the asymmetric distribution of the data, we determined the relationships based on the Spearman rank correlation. We used SPSS Statistics version 21.0(Armonk, NY, USA) to perform the analyses [63].  We also related the environmental variables with macrophyte and phytobenthic diatoms communities. Five groups of environmental variable stressor gradients were defined: nitrogen (N variables); phosphorous (P variables); land use (CLC variables); nutrients (N and P variables); and nutrients and land use (N, P, and CLC variables). Prior to the analyses, the phytobenthic diatoms data were transformed using the ln(x + 1) function, whereas the macrophyte data were not transformed. For the nutrient variables, we used log(x + 1) transformation to correct the heteroscedasticity in the data. Canonical Correspondence Analyses (CCAs) direct ordination techniques were carried out to analyse the relations between different groups of stressor gradients and the macrophyte and phytobenthic diatoms communities. Altogether, 14 CCAs were performed, seven with macrophytes and seven with phytobenthic diatoms communities (Table 3). Besides the differences in the two aquatic organism groups, the CCAs differed also in selected environmental variable group(s) and the forward selection procedure. In ten CCAs, an automatic forward selection procedure was used and manual in four CCAs. When a manual forward selection procedure was used, in the nutrient group, the phosphorous variables were selected first and then the nitrogen variables, whereas, in the nutrient and land use stressor group, the land use variables were selected first and then the nutrient variables. In the first step, we examined the contribution of each environmental variable to the explained variance of the species data (marginal effect), while, in the second step, we examined the remaining contribution of each variable added to the model when the other variables were already taken into account (conditional effects). Significant variables (p < 0.05) were selected using the Monte Carlo permutation test with 999 unrestricted permutations. These analyses were performed using the Canoco 5 software package [64].   We also related the environmental variables with macrophyte and phytobenthic diatoms communities. Five groups of environmental variable stressor gradients were defined: nitrogen (N variables); phosphorous (P variables); land use (CLC variables); nutrients (N and P variables); and nutrients and land use (N, P, and CLC variables). Prior to the analyses, the phytobenthic diatoms data were transformed using the ln(x + 1) function, whereas the macrophyte data were not transformed. For the nutrient variables, we used log(x + 1) transformation to correct the heteroscedasticity in the data. Canonical Correspondence Analyses (CCAs) direct ordination techniques were carried out to analyse the relations between different groups of stressor gradients and the macrophyte and phytobenthic diatoms communities. Altogether, 14 CCAs were performed, seven with macro-  In the last step, we built the composite stressor gradient. This was developed when running the Principal Component Analysis (PCA) in the Canoco 5 software package using the nutrient variables and CLC land cover variables [64].

Relationships between Nutrient Variables and Biotic Indices
The relationships of the River Macrophyte Index (RMI EQR) and the Trophic Index (TI EQR) with the nutrient variables varied significantly. The TI EQR was negatively weakly to moderately statistically significantly correlated (p < 0.05) with several nutrient variables (r = −0.33 to −0.59), while the RMI EQR was correlated (p < 0.05) with the summer mean total nitrogen values ( Table 4). Table 4. The relationship (Spearman rank correlation coefficients) between the nutrient variables and the ecological quality ratio of the River Macrophyte Index (RMI EQR) and the Trophic Index (TI EQR); n = 70. Statistically significant relationships are marked **-p < 0.01 and *-p < 0.05. The relationships between the stressor gradients and biotic indices were affected by the length of the gradient. The strength of the relationship varied with the length of the gradient and the selected biotic index. For the variable orthophosphate-maximum value (PO4 max), we got a stronger relationship with the RMI index at 0.3 and at the 0.2-mg/L total phosphorous-annual mean value (TP mean) (Figure 4). The whole gradient provided much lower and nonsignificant values. When we examined the relationships with the Trophic Index (TI), the values were significant only at gradient lengths > 0.3 mg/L of the TP mean and highest at the whole gradient. For the variable ammonium-maximum value (NH4 max), the strength of the relationship also varied with the length of the gradient ( Figure 5). The RMI EQR revealed no statistically significant relationships, while the TI EQR showed statistically significant relationships at gradient lengths ≥ 0.3 mg/L of the TP mean.

Relationships between Nutrient, Land Use Variables, and Biotic Communities
We tried to explain the variability of macrophytic and phytobenthic community structure by different simple and composite gradients, and we found out that the explained variabilities were very similar for the macrophytes and phytobenthic diatoms (Tables 5 and 6 and Figure 6). When we examined the explanatory potential of a single variable, we found out that the most variability in the macrophyte community was explained by the variable NO 2 mean. A somewhat lower share was explained by both the nitrogen (TN max, TN mean, NO 3 max, NO 2 max, and NO 3 mean) and phosphorus (TP mean) variables and the CLC group variables (Intensive agriculture-catchment and Natural-catchment) ( Table 2). For the phytobenthic diatoms community, the variability was best explained by the NO 2 mean variable. The slightly lower share of variability was explained by the CLC variables (Intensive agriculture-catchment and Natural-catchment), nitrogen variables (TN mean, NO 3 mean), and phosphorus variables (TP mean) ( Table 6).
The relationships between the stressor gradients and biotic indices were affected by the length of the gradient. The strength of the relationship varied with the length of the gradient and the selected biotic index. For the variable orthophosphate-maximum value (PO4 max), we got a stronger relationship with the RMI index at 0.3 and at the 0.2-mg/L total phosphorous-annual mean value (TP mean) (Figure 4). The whole gradient provided much lower and nonsignificant values. When we examined the relationships with the Trophic Index (TI), the values were significant only at gradient lengths > 0.3 mg/L of the TP mean and highest at the whole gradient. For the variable ammonium-maximum value (NH4 max), the strength of the relationship also varied with the length of the gradient ( Figure 5). The RMI EQR revealed no statistically significant relationships, while the TI EQR showed statistically significant relationships at gradient lengths ≥ 0.3 mg/L of the TP mean.

Relationships between Nutrient, Land Use Variables, and Biotic Communities
We tried to explain the variability of macrophytic and phytobenthic community structure by different simple and composite gradients, and we found out that the ex- Figure 5. Changes in the correlation coefficient (R) between the ammonium-maximum values (NH4 max) and River Macrophyte Index (RMI EQR) and Trophic Index (TI EQR), respectively, as a function of the stressor gradient length expressed as the total phosphorous-annual mean value (TP mean). Table 5. Percentage of macrophyte assemblage variance explained by each environmental variable before forward selection (λ1) and after forward selection (λa) within each stressor gradient group. Values of λa are given only for statistically significant explanatory variables (p < 0.05). LandNut-land use and nutrients, NutLand-nutrients and land use, P&N-phosphorus and nitrogen, Nut-phosphorous and nitrogen, N-nitrogen, P-phosphorous, and Land-land use.

Composite Stressor Gradient
In the last step, we built the composite stressor gradient using a PCA. Along the first and second axes of the PCA plot, we obtained different land cover and nutrient gradients. However, the nutrient variables were poorly correlated with the first or second PCA axes (Spearman rank correlation < 0.5), which reflects the independent gradients in our data set. The first PCA axis represented mainly a gradient of intensive agriculture (mostly in the subcatchment), while the second axis represented a gradient of natural areas and urban areas (in opposite directions) in the subcatchment (Figure 7). The relationships between the RMI EQR values, TI EQR values, and the PCA1 and PCA2 varied and also depended on the length of the gradient. The RMI EQR values were statistically significantly correlated only with PCA1 up to 0.3 mg/L of the total phosphorous mean annual value, representing the intensive agriculture in the subcatchment. Stronger relationships were observed for shorter gradients of lower values and decreased as the gradients extended towards higher values (Figure 8). The TI EQR values did not show any statistically significant correlations with any PCA axes of any stressor gradient length. When we applied different stressor gradient groups, we explained the greatest share of variability in the macrophyte and phytobenthic communities by the LandNut stressor group (22-26%), the least by phosphorus variables (phytobenthic diatoms, 9.5%; macrophytes, 11.7%), or nitrogen variables (macrophytes, 11.7%) ( Figure 6). Even though the gradient analyses revealed that land use and nutrients explain some of the variability of macrophyte communities, considering the data of the selected 70 samples, we did not find statistically significant correlations between the River Macrophyte Index (RMI) and nutrients. The gradient analysis revealed that the land use and nutrients groups explained different portions of the variability of macrophytes.

Composite Stressor Gradient
In the last step, we built the composite stressor gradient using a PCA. Along the first and second axes of the PCA plot, we obtained different land cover and nutrient gradients. However, the nutrient variables were poorly correlated with the first or second PCA axes (Spearman rank correlation < 0.5), which reflects the independent gradients in our data set. The first PCA axis represented mainly a gradient of intensive agriculture (mostly in the subcatchment), while the second axis represented a gradient of natural areas and urban areas (in opposite directions) in the subcatchment (Figure 7). The relationships between the RMI EQR values, TI EQR values, and the PCA1 and PCA2 varied and also depended on the length of the gradient. The RMI EQR values were statistically significantly correlated only with PCA1 up to 0.3 mg/L of the total phosphorous mean annual value, representing the intensive agriculture in the subcatchment. Stronger relationships were observed for shorter gradients of lower values and decreased as the gradients extended towards higher values (Figure 8). The TI EQR values did not show any statistically significant correlations with any PCA axes of any stressor gradient length. Water 2021, 13, x FOR PEER REVIEW 15 of 21 Figure 7. PCA ordination diagram of the first two ordination axes with the land use and nutrient variables (n = 70). Environmental variables correlated <0.5 with PCA1 or PCA2 were omitted. Urb-urban, Nat-natural, IAg-intensive agriculture, Cat-catchment, and SCat-subcatchment. See Table 1 for the environmental variable codes.   Table 1 for the environmental variable codes. Urb-urban, Nat-natural, IAg-intensive agriculture, Cat-catchment, and SCat-subcatchment. See Table 1 for the environmental variable codes.

Discussion
Nutrients are crucial for river biocenoses, supporting primary production, on one hand, while, on the other hand, their high concentrations compromise species relations within the community [65] and, thus, ecological status of the rivers. The response of macrophyte communities to nutrients differed significantly from that of diatom communities and so did the indices based on these two groups of organisms. This is possibly a consequence of the sampling time, since the macrophytes data used in our study were obtained at the peak growing season. The majority of studies relate to nutrients with macrophyte abundance only; however, we took into account species presence and abundance. For example, the study of O'Hare [66] revealed the increase of macrophyte biomass with P availability as a soluble reactive phosphorus and total phosphorus in water. The study of D'Aiuto [67] showed a significant correlation between stream nutrient enrichment (total phosphorus, soluble reactive phosphorus, and nitrate nitrogen) and macrophyte production. On the contrary, with epiphytic algae that obtains nutrients from water, an important aspect defining the response of the majority of macrophyte species to nutrients is also root uptake from the sediment. Carr [26] showed that, in streams, the macrophyte abundance increased linearly with increasing the sediment P concentration in water, while the study of Mebane [68] revealed that both the nitrogen contents in sediment and water were positively correlated to a macrophyte biomass. Thus, the response of the macrophyte community largely reflects the trophic conditions during their entire growth cycle and on the long-term accumulation of nutrients in the sediment, as shown by Clarke [28]. The absence of significant relationships of the RMI with the majority of the nutrient parameters is also a consequence of specific traits of species and communities, namely species life cycle, their growth dynamics, and their modes of nutrient uptake [69][70][71]. In the case of perennials with different underground storage organs, the impact of nutrients can be even more delayed and observed only over a long period [72]. An additional reason for the lack of correlation between the RMI with total phosphorus in our study could be the consequence of the relatively high total phosphorus concentrations in Slovenian watercourses, since, in our dataset, only 11% of the samples contained concentrations below 0.05 mg/L. The research of Fabris [73] applying the macrophyte Trophic Index (TIM) as a measure of the macrophyte community status revealed significant changes at relatively low concentrations of total phosphorus (up to 0.05 mg/L).
For the TI, we determined weak-to-moderate statistically significant correlations for the majority of the studied nutrient parameters. The significance of the correlations for the TI was expected due to the expansion mode of phytobenthos growth [11,74,75] and their quick response to nutrients [39].
The strength of the relationships between the RMI and TI and nutrients strongly depended on the length of the stressor gradient. The findings based on shorter stressor gradients with lower values (up to 0.3 mg/L of the TP mean) differed significantly from those based on longer ones. These inconsistent responses of the RMI relationships to nutrients may be a consequence of macrophytes that can function as a sink and a source of nutrients [76,77]. It was shown that the bulk of nutrients can be released from macrophyte beds during decay processes after disturbances and at the end of the vegetation period [78]. Besides that, macrophytes can also affect the nutrient balance by trapping fine organic and inorganic particles, enhancing the mineralization of organic matter through oxidation of the sediments, and altering the local environment, thus enabling P release due to the reducing conditions [79]. For the TI, statistically significant relationships with an annual mean total phosphorous value were obtained at gradient lengths ≥ 0.3 mg/L of the TP mean. For the variable maximum value of ammonium, the strength of the relationship also varied with the length of the gradient in the case of the RMI, but it showed an increase with the length in the case of the TI. Although some studies indicate a strong correlation between the stressor gradient and biotic indices/communities [80], the weak correlations obtained in our study were not against our expectations [61]. Weak correlations between the groups of stressor gradients (e.g., eutrophication) or individual stressor gradient (e.g., TP stressor gradient) and biotic indices confirmed that it is difficult to distinguish among the effects of individual stressor gradients or a group of stressor gradients on biotic communities in a multiple-stressor ecosystem. Moreover, different responses of biotic communities to the same stressor gradients are also possible [81]. Even in a case when we primarily assess the individual stressor gradient, there are always other events or influences that could change the conditions. For instance, removing riparian vegetation, which is a hydromorphological element, affects the water temperature, nutrient cycling, and consequently, the nutrient uptake and trophic conditions for biotic communities [82,83]. In addition, different environmental factors may act at multiple scales [84].
The explained variability of macrophytic and phytobenthic assemblage structures by different stressor gradient groups was similar for macrophytes and phytobenthic diatoms. The greatest share of variability was explained by the land use and nutrients (LandNut) stressor group and the least by the phosphorus variables and nitrogen variables. Thus, the unexplained portion of variability was still prevalent. Different studies revealed that, even with indices that focus on individual stressors (e.g., trophic status), authors rarely explain >50% of the variability of an index [85]. However, in all these cases, one does not evaluate only the trophic status but also the general degradation of the riverine ecosystem. To comprise these complex effects, we built a composite gradient using a PCA. Along the first and second axes of the PCA plot, we obtained different land cover gradients. The relationship with the land cover variables was expected, since many studies showed that the land cover presented a key factor affecting the aquatic community [15,22,86]. The first axis represented mainly a gradient of intensive agriculture, while the second axis a gradient of natural areas and urban areas. The gradient vectors of natural areas and urban areas pointed in opposite directions. The relationships between the RMI EQR values, TI EQR values, and the PCA1 and PCA2 varied and were dependent on the gradient lengths. In the case of the RMI EQR values, stronger relationships were observed for shorter gradients of lower stressor values and decreased as the gradient extended towards higher values. We obtained a significant relationship with PCA1 at low mean annual values of the total phosphorus, which was related to intensive agriculture in the subcatchment. With the TI EQR values, no significant relationships were determined. This may be due to the long-term response of macrophytes to their habitat conditions, including those related to land use changes, which may result in continuous sediment accumulation and nutrient input [87], while phytobenthos mainly responded to short-term changes of the nutrients in the water [1,8].
The observed similarities in the explained variability of the macrophytic and phytobenthic assemblage structures by different stressor gradient groups (nutrients and land use) and differences in the indices performances might indicate the individual assessment index limitations. Each index was developed with a specific dataset and tested against selected stressor gradient data. Usually, the stressor gradient applied in studies is for necessity, a simplification of the true stressor gradient to which biota are exposed. Despite all the limitations, biological assemblages offer more stable information about the ecological status of a given ecosystem than abiotic variables, although we cannot always perfectly relate them. Macrophytes and phytobenthos assemblages reflect varied faces of the same river ecosystems.

Conclusions
The responses of phytobenthic diatoms and macrophyte communities to stressors in an aquatic environment differed significantly and so did the responses of the biotic indices that were developed based on these two groups. The relationships between the nutrient variables on the RMI and TI revealed the effect of the length of the stressor gradient, namely a stronger relationship between the RMI and annual mean total phosphorous was found at smaller gradient lengths in comparison to the TI. With different stressor gradient groups, we explained similar but relatively low shares of macrophyte and phytobenthic diatom community variabilities. The greatest share of variability in the macrophyte and phytobenthic diatom communities was explained by the combination of land use variables and nutrients (LandNut stressor group) and the least by phosphorus variables and nitrogen variables. The relationships of the RMI and TI to the composite stressor gradient, which comprised multiple variables shaping the phytobenthic diatoms and macrophyte communities, varied and were dependent on the gradient lengths. The RMI EQR values showed a significant relationship with PCA1 at low mean annual values of the total phosphorus, while, with the TI EQR values, we obtained no significant relationships. The used assessment methods showed that macrophytes respond better to shorter stressor gradients of lower values, whereas phytobenthic diatoms respond better to longer stressor gradients and higher values.