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

Abstract

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.

1. 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 near-natural 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.

2. Materials and Methods

2.1. Study Area

The study was conducted in an area of two neighbouring countries: Slovenia with a total area of 20,273 km² and 4573 km of rivers with catchments larger than 10 km², and Croatia with a total area of 56,594 km² and 12,884 km of rivers with catchments larger than 10 km² (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² (80.8%) of the Slovenian territory and 35,101 km² (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. Main characteristics of sampled rivers with the number of sampling sites and samples.

Country	River	Eco-Hydromorphological Type	Catchment Size Range (km2)	Altitude Range (m a.s.l.)	No. Sites (Samples)
Slovenia	Drava	Intermountain	11,720–13,091	253–338	6 (14)
Slovenia	Drava	Lowland-braided	13,189–15,079	178–236	12 (23)
Slovenia	Mura	Lowland-braided	9784–10,506	165–246	11 (21)
Slovenia	Sava	Intermountain	4946–5203	191–222	3 (7)
Slovenia	Sava	Lowland-deep	7151–7655	154–191	4 (8)
Slovenia	Sava	Lowland-braided	7782–10,411	132–139	3 (9)
Croatia	Mura	Lowland-braided	10,930–10,930	153–153	1 (1)
Croatia	Mura	Lowland-deep	11,731–11,731	141–141	1 (1)
Croatia	Drava	Lowland-braided	14,363–31,038	122–190	2 (3)
Croatia	Drava	Lowland-deep	33,916–39,982	81–100	4 (4)
Croatia	Sava	Lowland-braided	10,997–12,316	113–132	2 (2)
Croatia	Sava	Lowland-deep	12,884–64,073	74–91	6 (7)
Croatia	Kupa	Lowland-deep	9184–9184	92–92	1 (2)
Croatia	Una	Lowland-deep	9368–9368	94–94	1 (2)
Total			4946–64,073	74–338	57 (104)

). 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 gravel-bed 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² and altitudes between 74 and 338 m were included in this analysis.

2.2. 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 land-use, and water quality (

Table 2. Groups of environmental variables with their Median, Min (minimum) and Max (maximum) values.

Environmental Variable	Unit	Variable Group	Code	Median (Min-Max)	Transformation
Catchment size	km2	Typology	C_size	11,720.1 (4945.8–64,073)	log(x + 1)
Depth mean	Classified 1–2	Typology	Depth	1 (1–2)
Slope	(%)	Typology	slope	0.9 (0–3.6)
Altitude	(m a.s.l.)	Typology	altitude	190 (74–338)
Substrate	Classified 1–3	Typology	substratum	3 (1–3)
Water temperature	°C	Water quality	T	12 (8.6–20.5)	log(x + 1)
pH		Water quality	pH	8 (7.7–8.3)	log(x + 1)
Conductivity	µS/cm	Water quality	cond	339.5 (262–517)	log(x + 1)
Oxygen concentration	mg O2/L	Water quality	DO	9.5 (7.4–11.5)	log(x + 1)
Oxygen saturation	(%)	Water quality	DOsat	89 (79.5–104.2)	log(x + 1)
Total suspended solids	mg/L	Water quality	TSS	7 (2.4–56)	log(x + 1)
Chemical oxygen demand (K2Cr2O7)	mg O2/L	Water quality	COD	5.7 (2.5–19.6)	log(x + 1)
Biochemical oxygen demand (5 days)	mg O2/L	Water quality	BOD5	1.2 (0.7–2.9)	log(x + 1)
Orthophosphate	mg P/L	Water quality	PO4	0 (0–0.3)	log(x + 1)
Total nitrogen	mg N/L	Water quality	Ntot	1.6 (0.7–2.6)	log(x + 1)
Ammonium	mg N/L	Water quality	NH4	0 (0–0.3)	log(x + 1)
Nitrite	mg N/L	Water quality	NO2	0 (0–0.1)	log(x + 1)
Nitrate	mg N/L	Water quality	NO3	1.4 (0.5–2)	log(x + 1)
Urban land use	(%)	Land use	C_urb	3.1 (0.9–5)	arcsin(sqrt x)
Natural and semi-natural land use	(%)	Land use	C_nat	70.8 (55–78.6)	arcsin(sqrt x)
Non-intensive agriculture land use	(%)	Land use	C_agrE	12.1 (10.1–24.9)	arcsin(sqrt x)
Intensive agriculture-tilled land use	(%)	Land use	C_agrI1	3.6 (0.7–14.7)	arcsin(sqrt x)
Intensive agriculture-non-tilled land use	(%)	Land use	C_agrI2	7.2 (4.5–20.9)	arcsin(sqrt x)
Discharge	m3/s	Hydromorphology	Q	119.2 (5.2–824)	log(x + 1)
Mean annual discharge	m3/s	Hydromorphology	NQ	114.7 (6.6–648)	log(x + 1)
Lowest annual discharge	m3/s	Hydromorphology	MQ	216.7 (9–998)	log(x + 1)
Highest annual discharge	m3/s	Hydromorphology	HQ	483.3 (31.1–1530)	log(x + 1)
River habitat quality index	Total score 1	Hydromorphology	RHQ	218.6 (43.5–324.3)
River habitat modification index	Total score 1	Hydromorphology	RHM	28.6 (0–116)
Hydrological modification index	Total score 1	Hydromorphology	HLM	0.8 (0–1)
Hydromorphological modification index	Total score 1	Hydromorphology	HMM	0.7 (0–1)
Hydromorphological quality and modification index	Total score 1	Hydromorphology	HQM	0.7 (0–1)

¹ Calculated score of individual features according to the SIHM method [20].

). 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(K₂Cr₂O₇), BOD₅, 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] using ArcGIS version 10.2.1 (Esri Corp., Redlands, CA, USA). The categories were combined into five land use variables: urban land use (CLC class 1), natural and semi-natural land use (CLC classes 3, 4, 5), non-intensive agriculture land use (CLC categories 2.3.1, 2.4.3, 2.4.4), intensive agriculture land use on arable land (CLC categories 2.1), and intensive agriculture land use on non-arable land (CLC categories 2.2, 2.4.1, 2.4.2). Nine hydromorphological (HM) variables were selected; four discharge parameters and five HM indices of the Slovenian hydromorphological (SIHM) assessment method [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 eco-hydromorphology 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.

2.3. 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² 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).

2.4. 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 (R_Sp) 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.

3. Results

3.1. Relationships between the Variables

Spearman rank correlation (R_Sp) resulted in several statistically significant relationships (P < 0.05) between pairs of environmental variables (Appendix C, Appendix D, Appendix E and Appendix F). Strong correlations (|R_Sp| > 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 < |R_Sp| < 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 (|R_Sp| < 0.50) or insignificant.

3.2. 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. Percentage of benthic invertebrate assemblages variance explained by each environmental variable’s independent (before forward selection) and conditional effects (after forward selection within each environmental variable group).

		Before FS	After FS Groups
Environmental Variable	Variable Group	ƛ	ƛ	P	F
Depth	Typology	0.16	0.16	0.001	2.86
Altitude	Typology	0.14	0.14	0.001	2.51
Slope	Typology	0.11	0.07	0.122	1.3
Catchment size	Typology	0.1	0.09	0.001	1.73
Substratum	Typology	0.07	0.09	0.013	1.63
Conductivity	Water quality	0.22	0.22	0.001	3.93
Nitrogen-total	Water quality	0.14	0.06	0.444	1.01
Nitrate	Water quality	0.13	0.05	0.358	1.05
Nitrite	Water quality	0.13	0.06	0.273	1.09
Orthophosphate	Water quality	0.13	0.09	0.001	1.75
Ammonia	Water quality	0.11	0.09	0.011	1.68
COD	Water quality	0.1	0.07	0.066	1.29
Temperature	Water quality	0.08	0.07	0.04	1.34
BOD5	Water quality	0.08	0.06	0.12	1.21
Dissolved oxygen saturation	Water quality	0.08	0.08	0.031	1.41
Dissolved oxygen concentration	Water quality	0.07	0.07	0.1	1.2
pH	Water quality	0.07	0.05	0.412	1.02
Total suspended solids	Water quality	0.06	0.07	0.033	1.36
Natural and semi-natural land use	Land use	0.12	0.12	<0.0001	2.19
Intensive agriculture-non-tilled land use	Land use	0.11	0.12	<0.0001	2.12
Intensive agriculture-tilled land use	Land use	0.11	0.06	0.187	1.16
Non-intensive agriculture land use	Land use	0.11	0.09	0.005	1.72
Urban land use	Land use	0.08	0.1	0.002	1.81
Hydrological modification index	Hydromorphology	0.18	0.18	0.001	3.29
Hydromorphological quality and modification index	Hydromorphology	0.17	0.08	0.011	1.62
River habitat quality index	Hydromorphology	0.15	0.08	0.001	1.61
Hydromorphological modification index	Hydromorphology	0.15	0.08	0.015	1.45
Discharge	Hydromorphology	0.1	0.08	0.015	1.51
River habitat modification index	Hydromorphology	0.1	0.1	0.002	1.73
Highest annual discharge	Hydromorphology	0.08	0.07	0.039	1.29
Lowest annual discharge	Hydromorphology	0.08	0.08	0.015	1.51
Mean annual discharge	Hydromorphology	0.08	0.05	0.528	0.96

). 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.

3.3. 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.

3.4. 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.

4. 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 near-natural 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.

5. Conclusions

• 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.