Disentangling the Effects of Multiple Stressors on Large Rivers Using Benthic Invertebrates—A Study of Southeastern European Large Rivers with Implications for Management

: Predicting anthropogenic actions resulting in undesirable changes in aquatic systems is crucial for the development of effective and sustainable water management strategies. Due to the co ‐ occurrence of stressors and a lack of appropriate data, the effects on large rivers are difficult to elucidate. To overcome this problem, we developed a partial canonical correspondence analyses (pCCA) model using 292 benthic invertebrate taxa from 104 sites that incorporated the effects of three stressors groups: hydromorphology, land use, and water quality. The data covered an environmental gradient from near ‐ natural to heavily altered sites in five large rivers in Southeastern Europe. Prior to developing the multi ‐ stressor model, we assessed the importance of natural characteristics on individual stressor groups. Stressors proved to be the dominant factors in shaping benthic invertebrate assemblages. The pCCA among stressor ‐ groups showed that unique effects dominated over joint effects. Thus, benthic invertebrate assemblages were suitable for disentangling the specific effect of each of the three stressor groups. While the effects of hydromorphology were dominant, both water quality and land use effects were nearly equally important. Quantifying the specific effects of hydromorphological alterations, water quality, and land use will allow water managers to better understand how large rivers have changed and to better define expectations for ecosystem conditions in the future.


Introduction
It is recognised that large rivers are economically important, but they also provide various ecosystem services and hence require sustainable management. The European Water Framework Directive [1] requires the identification of significant anthropogenic pressures and the assessment of their impacts on water bodies. Thus, we need to correctly predict human activities that create unacceptable impacts on large rivers. While the sources of stress in large rivers are numerous [2,3], little is known about the prevalence, spatial patterns, interactions with the natural environment and co-occurrence of stressors and their effects [4]. The effects of multiple stressors are difficult to predict due to the complexity of the interactions among stressors [5,6]. Thus, the effects of the individual stressor may be masked by the presence of other stressors.
Human pressures and land use patterns have long been recognised as a threat to the functioning and ecological integrity of aquatic ecosystems, as impacts on habitats, water quality, and biota involve complex pathways, e.g., [7,8]. High amounts of pollutants and nutrients have been discharged into large rivers as a result of industrial development, urbanization and intensive agriculture [8]. During the 19th and 20th centuries, stream regulation transformed large rivers to allow for navigation and power generation at the expense of habitat loss [9,10]. Large rivers became impounded, and their channels straightened and separated from oxbow lakes by levees to protect human settlements against floods [11][12][13]. These activities generally reduced longitudinal connectivity and connectivity between the main channel and adjacent floodplain channels [14], disturbing the natural gradients of chemical and physical parameters along large river courses were disturbed.
Aquatic communities are altered on a relatively predictable gradient from natural, e.g., undisturbed or minimally disturbed conditions, to severely altered conditions [15]. Ecological studies of large rivers are usually limited to individual rivers and are rarely based on data along the whole environmental gradient from near-natural to heavily altered sites. The reason might be that in individual large rivers, especially of developed countries including Europe, few near-natural remain [16]. However, in Southeastern Europe, despite large, altered stretches, the large rivers contain some of the last natural, free-flowing stretches in Europe. The Kupa and Una Rivers are in near natural conditions along their entire courses. The middle and lower stretches of the Drava and Mura Rivers are rare examples of unregulated, very large European rivers. Major lower sections of the Sava River still exhibit a relatively natural geomorphic structure and hydrological regime and are fringed by large protected wetlands.
Certain natural characteristics (e.g., catchment characteristics, depth, channel pattern) play a role in structuring benthic invertebrate communities, even in large river and at the regional scale e.g., [17]. Thus, the differences in these characteristics across large rivers must be accounted for before the impacts of stressors can be examined. Aside from natural conditions, hydromorphological alterations (the concept of 'hydromorphology' is a term introduced by the EC Water Framework Directive [1] that includes hydrological, morphological, and river continuity characteristics), land use, and water quality profoundly affect benthic invertebrates in rivers. Understanding the specific and joint effects of these stressors is of critical importance for developing effective river basin management plans to shape environmental policy.
In this study, we examined the unique and joint effects (two or more factors) of natural factors and major stressors (hydromorphology, land use, and water quality) on the invertebrate fauna of Southeastern European large rivers using the data along the entire environmental gradient from nearnatural sites up to heavily altered sites. The term stressor(s) refers to variable(s) of anthropogenic landscape changes and local abiotic stream conditions that reflect human activities, and herein is used in this sense. Natural factors not influenced by anthropogenic disturbance are referred to using the term typology. We posed three general hypotheses regarding benthic invertebrate responses to natural factors (typology) and major stressors: (1) Stressors and natural factors play a key role in structuring benthic invertebrate communities in the large rivers of a certain region (e.g., Southeastern Europe), thus differences in natural characteristics must be accounted for before the impacts of stressors can be isolated.
(2) Hydromorphology, land use, and water quality have distinct individual effects on structuring stream benthic macroinvertebrate assemblages.
(3) Specific stressor effects of hydromorphology, land use, and water quality are more important than their joint effects in structuring the benthic macroinvertebrate assemblages of large rivers and thus benthic invertebrates can be used to disentangle the effects of these stressors on large rivers.

Study Area
The study was conducted in an area of two neighbouring countries: Slovenia with a total area of 20,273 km 2 and 4573 km of rivers with catchments larger than 10 km 2 , and Croatia with a total area of 56,594 km 2 and 12,884 km of rivers with catchments larger than 10 km 2 (Figure 1). The rivers in each of the countries belong either to the Danube or Adriatic River Basin, though this study included only rivers of the Danube River Basin. The Danube River Basin covers 16,381 km 2 (80.8%) of the Slovenian territory and 35,101 km 2 (62%) of Croatia. The landscape within this basin is diverse in altitude and slope and features different river section types [18]. This study was limited to five major rivers of the Danube River Basin: the Drava with Mura, and the Sava with its tributaries Kupa/Kolpa and Una (Table 1). The Sava and Drava Rivers are among the largest discharge tributaries of the Danube River (1st and 4th, respectively) and represent some of the best-preserved rivers in Europe in terms of their biological and landscape diversity. The Sava River springs in Slovenia as a gravelbed river under Alpine influences, the channel in Slovenia changes from simple straight to braided, before gaining its meandering course downstream of Zagreb and continuing to its mouth in Belgrade (Serbia). The Sava River is considered by nature conservationists and scientists to be one of the crown jewels of European nature [19]. The Drava River crosses ecoregions from high Alpine mountains to the Pannonian-Illyrian plain and features all typical fluvio-morphological river types from straight to braided to meandering channels. The lower Drava with the lower Mura River constitutes a 380 km free-flowing and semi-natural watercourse and represents one of the last remaining continuous, riverine landscapes in Central Europe [18]. Only stretches with a catchment area from between 5000 and 64,000 km 2 and altitudes between 74 and 338 m were included in this analysis.

Environmental Variables
The sampling sites cover near-natural to highly disturbed conditions, reflecting the various disturbance levels caused by different stressors, e.g., hydromorphological alteration, catchment landuse, and water quality ( Table 2). A total of 34 environmental variables were measured or calculated and classified into four groups: typology (natural), hydromorphology, land use, and water quality. The data for the five typology variables were obtained from the GIS database, the hydrological databases of the Slovenian Environment Agency (ARSO) and Croatian Meteorological and Hydrological Service (DHMZ), and from field analyses. Altitude and slope were calculated using a digital elevation model with 5 m accuracy. The natural predominant substrate was classified into three classes, as the smaller fraction (psammal-1), small to medium fraction (psammal/akal-2) and larger fraction (lithal-3), and the mean depth at low water level was defined as 1 when <1.5 m and 2 when >1.5 m.
Physical and chemical data were obtained monthly or at least four times a year (each season) from the national surface water monitoring programmes. In these analyses, only those 13 parameters were considered where data were available for all selected sites (Table 2): conductivity, pH, oxygen concentration, oxygen saturation, water temperature, COD(K2Cr2O7), BOD5, orthophosphate, total nitrogen, ammonium, nitrite, nitrate, and total suspended solids. In the analyses, the median of data gathered for each parameter in the year of benthic invertebrate sampling was used.
Land use variables were defined from the share of land use categories at the catchment scale, extracted from Corine Land Cover (CLC) data [20] [17,21,22]. The SIHM method was applied to examine habitat quality, habitat modifications, and the influence of main upstream barriers/impoundments were considered. First, a river habitat survey [23,24] was performed once for each sampling site and the data was used to calculate the two morphological indices [21,22]: river habitat quality index (RHQ), and river habitat modification index (RHM). Normalised values (converted to a common scale of 0-1; RHQnor, RHMnor; [17]) were used. We first defined the ecohydromorphology types of the considered river stretches according to [17] (Table 1, Figure 1). The RHM index was normalised using the same values for all river types; a reference value and a lower anchor of 0 and 112, respectively. RHQ values were normalised using type specific reference values. For intermountain and lowland-deep eco-hydromorphology river types, a reference value of RHQ = 237 was used, whereas for lowland-braided the RHQ was set at 327. The lower anchor was the same (RHQ = 116) for all river types. The data on impoundments recorded in the catchment of each sampling site was used to calculate the hydrological modification index (HLM; [17,21]). Combining the indices RHQnor, RHMnor and HLM two HM indices were calculated: hydromorphological modification index (HMM) and hydromorphological quality and modification index (HQM) [17,21,22]. Hydrological variables were obtained from the available data on discharge from national monitoring gauging stations (ARSO, DHMZ). In addition to the mean daily value of discharge measured on the day of benthic invertebrate sampling (Q), the mean annual discharge (MQ), the lowest annual discharge (daily average; NQ), and the highest annual discharge (daily average; HQ) were also calculated for the sampling-year period.

Benthic Invertebrates
Biological data were obtained as part of the WFD monitoring and assessment system development programmes in Slovenia and Croatia between 2005 and 2011. In total, 104 samples were collected at 57 sites: 39 sites (82 samples) in Slovenia and 18 sites (22 samples) in Croatia (Figure 1, Appendix A). Some sites were sampled several times, but not more than once per year. Benthic invertebrates were collected during low to medium discharge using a multi-habitat sampling approach. Samples were collected in the wadeable part (up to 1.2 m) of the main channel or in the littoral zone of the impoundments to a depth of 1 m using a hand net (frame 25 × 25 cm, mesh-size: 500 μm). On each occasion, at every site, 20 sub-sampling units with a total sampling area of 1.25 m 2 were taken along a 100-250 m river stretch. The sampling procedure in Slovenia followed the standardized Slovenian river bioassessment protocol [17,25,26]. Twenty sampling units were selected in proportion to the coverage of the microhabitat types [17,24]. Microhabitat types were defined as the combination of substrate and flow type with at least 5% coverage. The channel substrate of each sampling site was classified according to [27], and flow characteristics according to [27,28]. Sampling units were pooled, preserved with 96% ethanol in the field and transferred to the lab for further processing. Each sample was sub-sampled, and the benthic organisms from a quarter of the whole field sample were identified and enumerated [29]. In Croatia, samples were collected according to the AQEM sampling strategy [27]. A total of 20 sampling units were sampled from representative substrates (i.e., substrates >5% coverage in the sample reach). At sampling sites with homogenous substratum (sand and other soft sediments) 10 sub-sampling units were taken instead of 20 (five sampling sites). In such cases, the sample was taken by pushing the hand net through the upper part (2-5 cm) of the substratum. The sampling units were pooled, preserved with 96% ethanol in the field, and transferred to the lab for further processing. In 2006, a more elaborate sub-sampling design was used, and habitat (substrate)-specific subsampling units were pooled and analysed as separate samples. In the lab, at least 1/6 of the sample was sorted until the minimum targeted number of 500 (habitat-specific samples) or 700 individuals (multi-habitat samples) was reached. Benthic invertebrates were identified usually to the species and genus level, though Oligochaeta and Diptera were identified to the (sub) family and genus level (Appendix B).

Data Analyses
Direct ordination techniques were carried out to analyse associations among environmental variables and between different groups of environmental variables and benthic invertebrate assemblages. These analyses were performed using Canoco 5 [30]. Benthic invertebrate data were transformed (ln(x + 1)) prior to analysis. In addition, some environmental variables were transformed prior to the analyses to approximate the normal distribution [31] (Table 2). Catchment size, water quality variables, and hydrological variables were transformed using log(x + 1), whereas land use data (proportional data) were transformed using arcsin(sqrt x). Spearman rank correlation coefficients (RSp) were calculated between all pairs of environmental variables using SPSS Statistics version 21.0 [32]. The rationale was to identify associations among the analysed groups of variables and to compare them among different datasets. Since sampling season of benthic invertebrates differed among the samples, prior to performing the direct ordination analysis, the importance of the temporal variable represented by the sampling day in a year was tested. As the temporal variable explained only a low percentage in the variance of the benthic invertebrate dataset, in comparison to the environmental variables, it was not included in the further analyses.
To determine the compositional gradient length the invertebrate data were analysed using detrended correspondence analysis (DCA; [33]). Since these gradient lengths were greater than two standard deviations, we assumed unimodal species responses and, thus, canonical correspondence analysis (CCA; [34]) and partial canonical correspondence analysis were applied [35]. For the first overview of the relationship between the environmental variables and benthic invertebrate data, a CCA analysis with an automatic forward selection routine was applied to all environmental variables. This process specified the effects that each environmental variable added to the explained variance of the species data (marginal effects) and the remaining effect that each variable added to the model once when other variables had already been loaded (conditional effects) [34]. Significant variables were selected with forward selection routine, using the Monte Carlo permutation test with 999 unrestricted permutations. The same procedure was then applied within variable groups (i.e., typology, hydromorphology, water quality, land use). The selected variables were used for partitioning the explained variance among benthic invertebrate assemblages using partial CCA (pCCA). This test allows for the investigation of the effects of one variable group, while eliminating the effects of other variable groups, and hence the partitioning of the variance into unique and joined effects of variable groups. The total explained variance among benthic invertebrate assemblages with forward selected environmental variables from three groups was partitioned into (i) the variance uniquely explained by each variable group, (ii) the variance explained by combined effects of each pair of variable groups, and (iii) the variance explained by combined effects of all three variable groups together.

Relationships between the Variables
Spearman rank correlation (RSp) resulted in several statistically significant relationships (P < 0.05) between pairs of environmental variables (Appendices C-F). Strong correlations (|RSp| > 0.70) among the variables of different stressor-groups were rare; natural and semi-natural land use related positively to altitude and negatively to conductivity, whereas alternatively non-intensive agriculture land use and intensive agriculture-non-tilled land use were negatively correlated with altitude and positively with conductivity. Other strong correlations were observed within all the stressor-groups, with the exception of the typology group. In the hydromorphology group, strong positive correlations were observed among hydrological variables and among indices (HLM, HMM, HQM). Several variable pairs of different stressor-groups showed moderate correlations (0.50 < |RSp| < 0.70). The lowest number of moderate and strong correlations was observed between the groups hydromorphology and water quality or land use. However, the indices HLM, HMM, and HQM showed a moderate positive correlation with conductivity and a negative correlation with natural and semi-natural land use. In the typology group, hydrological indices showed a moderate positive correlation with mean depth and catchment size and negative with slope. Most of the pairwise correlations were weak (|RSp| < 0.50) or insignificant.

Benthic Invertebrate Response to Environmental Variables
The total amount of variance (inertia) in the species data was 5.808, including the 104 sites and 292 benthic invertebrate taxa (Appendix B). The total explained variance in the dataset, including all 32 environmental variables, was 2.366 (41%). When tested individually, the highest explanatory power was observed for conductivity (0.22) (Table 3). Additionally, each of the 11 other variables showed more than 50% explanation power (>0.11) of the best explanatory variable. The variables of all four groups showed considerable explanation power. The hydromorphology group was represented with four indices (HLM index, HDM index, RHQ index and HMM index). An additional four eutrophication variables represented the water quality group (nitrogen-total, nitrate, nitrite, orthophosphate). Typology was represented by depth and altitude, whereas land use by natural and semi-natural land use. The remaining other 21 variables exhibited a weaker explanation power.
Testing each explanatory group individually, 22 of 32 environmental variables significantly contributed to the explained variance (Table 3). Each variable group comprised four to eight forward selected variables, used in the variance partitioning. The highest number of selected variables was in the hydromorphology group (eight out of nine), followed by the water quality group (six out of 13), typology group (four out of five) and land use group (four out of five). In the hydromorphology group, a combination of alteration indices and hydrological conditions was observed. For the water quality group, the selected variables reflected eutrophication, organic pollution and some other human activities. In the typology group, a combination of catchment conditions (catchment size, altitude) and instream conditions (depth, substrate) was observed. The land use group reflected urbanisation, agriculture, and other non-natural land use.

Variance Partitioning between Typology and Stressor-Groups
Variance partitioning between the typology group and an individual stressor-group revealed the unique effects of the stressor-groups, explaining from 36% (land use) to 53% (hydromorphology) of the benthic invertebrate assemblages explained variability (Figure 2). The joint effects (% of the total explained variance) of each stressor group and the typology group were relatively small (8-20%) in comparison to pure stressor effects. Joint effects always represented <30% of the stressor-group total explained variability.

Variance Partitioning of Three Stressor Variable Groups of Environmental Variables
Variance partitioning was run with 18 variables, after a forward selection routine for each variable group separately. Clearly, the unique effects of variable groups were more important in explaining the variation in the benthic invertebrate composition than joint effects (84% and 16% of the explained variance, respectively, Figure 3). The highest share (36%) was explained by the hydromorphology group, followed by water quality (27%) and land use group (21%). The explanatory power of any joint effect was much smaller where interaction between water quality and hydromorphology groups was most important, accounting for 10%. Other interactions between group pairs were less important (≤4%) and the joint effects of all three variable groups explained only 2% of the variation.

Discussion
Centuries of human activities including water pollution and habitat alterations have profoundly altered most large rivers and their aquatic assemblages. Conversion of native forests to agricultural and urban uses has increased concentrations of pollutants (e.g., nutrients), as well as habitat changes [36]. Only a few large European rivers still have stretches that appear to remain in their natural conditions. This study examined some of these stretches in Southeastern Europe, including natural and degraded large river stretches enabled us to cover the whole environmental gradient from nearnatural up to heavily altered sites. Thus, our sampling sites exhibit a wide range of chemical and physical factors reflecting differences in habitat characteristics, land use and water quality.
Large rivers are unique ecosystems. Although they share some abiotic and biotic commonalities, certain natural characteristics play a role in structuring benthic invertebrate communities of large rivers, even within the same region (e.g., [17]). Bonada et al. [37] stated that isolating the natural variability along a large river course from the influence of water pollution, land use, and hydromorphological (HM) alterations is difficult due to their confounding effects. This study confirmed the presence of certain joint effects of natural characteristics and individual stressors, though the joint effects were found to be less conspicuous than the specific stressor effects. Our results indicated that regional data of large rivers can be pooled and stressor effects isolated, partly supporting the hypothesis that differences in natural characteristics must be accounted before the impacts of stressors can be determined. We found that, for large rivers, the joint effects also depend on the stressor group. Water quality showed the highest joint effects with typology reflecting that water pollution impacts depend on the natural characteristics of the large river. For example, the effects of nutrients are more evident in large rivers with slower water flow and higher water temperature [38,39]. Hydromorphological alterations showed the lowest joint effects with typology, which might reflect that HM alterations similarly change benthic habitats and communities of large rivers in way in all large rivers. This supports the findings of minimal differences in responses of benthic invertebrate assemblages to HM alterations among large river types [17]. Land use showed an intermediate joint effect, possibly reflecting the combination of effects on water quality (e.g., eutrophication) and habitat characteristics (e.g., sedimentation) [25,36]. Due to the presence of multiple stressors and the lack of appropriate data, effects on large rivers can be difficult to elucidate [17,40]. In this study, pure stressor-specific effects on benthic invertebrates were heavily dominant over the joint effects, and hence, benthic invertebrate assemblages be useful in disentangling the effects of hydromorphology, water quality, and land use. Nevertheless, substantial joint effects were observed between HM alterations and water quality. The combination of both these stressors likely exerts substantial change in benthic invertebrate communities in large rivers. It is often observed that HM alterations (e.g., water abstraction, damming) also lead to water quality issues (e.g., eutrophication) [39,41,42]. The joint effects of the other two stressor groups and all three stressor groups together were small. It is known that land use changes impact water quality and HM conditions, thus, substantial joint effects could be expected. However, it seems that in Southeastern Europe, land use is not so intensive to severely influence water quality and/or HM conditions of large rivers. Moreover, river damming and channelling are key HM pressures impacting large river benthic habitats in the region [43,44]. We found that the effects of HM alterations were more significant than those of water quality and land use changes. Although the sequence and timing of individual stressor effects were not determined in this study, there is evidence that water pollution was most important during the early to mid-20th century [45,46]. Changes in land use have long been present but are intensifying, whereas HM alterations have become more significant in recent decades. Physical pressures have been identified as major causes for a potential failure of water bodies to meet the Water Framework Directive environmental objectives [4,47]. Nevertheless, this study showed that special attention should also be given to the effects of water pollution and land use.
Modelling the relationship between biological communities and environmental parameters has played an increasingly important role in ecology [48]. Such a predictive approach can lead to a better understanding of how the species composition can potentially be affected by human pressures, and is especially promising for use in conservation planning and resource management [49]. Partial canonical correspondence analyses (pCCA), proved to be a useful method for disentangling the effects of addressed stressors. However, there are limitations since the results also depend on selected variables, length of the gradient, and correlation among variables (e.g., [25,50,51]). In this study, the stressor variables were selected according to the reported impacts on large river aquatic communities [8,17,41]. In the water quality group, several other parameters could have been selected, though we chose the most relevant parameters from the eutrophication and organic pollution group. Strong correlations were found between environmental variables, but only within the stressor group. Therefore, it was possible to isolate the effects of different stressors on benthic invertebrates. The long environmental gradient of sites from near-natural conditions up to heavily altered sites is crucial for building a reliable model. It is also important to view the large river community dynamics not only in the context of environmental variables, but also in biotic interactions [52]. Alien species in particular might influence benthic invertebrate community responses (e.g., [50]), though this was not an issue in this study, as recorded alien species (e.g., Corbicula fluminea, Dreissena polymorpha, Dikerogammarus villosus, D. haemobaphes, Jaera istri) usually represented less than 5% of the benthic invertebrate assemblages' sample composition.
Understanding the impact of water pollution, hydromorphology, and land use change on the ecological status and ecosystem services is essential for developing effective river basin management plans (RBMPs) and shaping future environmental policy. Setting appropriate measures will enable environmental objectives to be achieved (e.g., good ecological status according to Water Framework Directive [1]). Relationships have previously been defined between the biota and water quality [53][54][55] what resulted in active river management for water quality improvement [56]. We showed that water quality issues still exist in large rivers and their effects also interact with HM alterations and land use. HM alterations are the dominant stressor in rivers throughout Europe [4,47], and many studies consider only HM alterations (e.g., [17,22,[57][58][59][60]). We showed that in addition to HM alterations and water quality, land use impacts on benthic invertebrates are substantial. For example, increased urbanisation and intensive agriculture severely impact benthic invertebrate assemblages and the integrity of large rivers. Therefore, all major stressors need to be addressed and their effects disentangled to ensure implementation of sustainable river basin management strategies (e.g., Water Framework Directive). The integration of environmental objectives in sectoral policies (e.g., Common Agricultural Policy, Floods Directive, renewable energy, Natura 2000) having direct or indirect impacts on rivers and their catchments might help to achieve environmental objectives (e.g., good ecological status). Since most large European rivers have catchments that cross international borders, cooperation among countries is critical in planning and implementing management strategies. However, differences in development level, public opinion, and historical and political constraints can hinder attempts to achieve these common environmental objectives. Southeastern Europe is facing a range of development challenges, including the planning of new hydroelectric power plants, ongoing intensive urbanisation, and intensifying agriculture [13,25,44].
Public understanding of the importance of water quality, habitat conditions and land use in structuring aquatic assemblages in large rivers could provide a basis for greater support of effective large river protection and sustainable management efforts. However, the management agencies of Southeastern Europe need to change their paradigm of river water quality to the ecological quality of the river ecosystem, thereby supporting activities that would prevent large river deterioration as was observed in many parts of the world.


We disentangled the specific effects of hydromorphology, water quality, and land use using benthic invertebrate assemblages.  Joint effects of stressors and natural factors on benthic invertebrate assemblages depend on the stressor group.  Stressors proved to be the dominant factors in shaping benthic invertebrate assemblages of Southeastern Europe large rivers. Effects of hydromorphology dominated over water quality and land use effects, though these were still substantial. Thus, all major stressors need to be addressed and their effects determined for the implementation of the sustainable river basin management strategies.  Management agencies in Southeastern Europe need to change their paradigm from river water quality to the ecological quality of the river ecosystem, thereby supporting activities that will prevent large river deterioration. Funding: This research was conducted within the Slovenian-Croatian bilateral project "Benthic invertebrate based ecological status assessment of large rivers with management goals focused on hydromorphological alterations" (to G.U. and Z.M). The study was also supported by the Ministry of the Environment and Spatial Planning of the Republic of Slovenia as a part of the National Programme for the Implementation of the EU Water Framework Directive, conducted at the Institute for Water of the Republic of Slovenia.
Acknowledgments: Authors gratefully acknowledge the project team members for their assistance both in the field and in the laboratory. We thank the Croatian Meteorological and Hydrological Service and Croatian Waters and Slovenian Environment Agency for providing hydrological and water physico-chemical data. We thank Maja Kerovec for assisting with GIS analysis.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Appendix C Table A3. Statistically significant Spearman's correlation coefficients (RSp) for the combinations of typology variables with variables from all groups (**P < 0.01, *P < 0.05); |RSp| > 0.50 are in bold. See Table 2 for environmental variable codes. Med-median.  Table A4. Statistically significant Spearman's correlation coefficients (RSp) for the combinations of water quality variables with variables from all groups (**P < 0.01, *P < 0.05); |RSp| > 0.50 are in bold. See Table 2 for environmental variable codes. Med-median.  Table A5. Statistically significant Spearman's correlation coefficients (RSp) for the combinations of land use variables with variables from all groups (**P < 0.01, *P < 0.05); |RSp| > 0.50 are in bold. See Table 2 for environmental variable codes. Med-median.  Table A6. Statistically significant Spearman's correlation coefficients (RSp) for the combinations of hydromorphology variables with variables from all groups (**P < 0.01, *P < 0.05); |RSp| > 0.50 are in bold. See Table 2 for environmental variable codes. Med-median.