Taxonomic and Functional Responses of Benthic Macroinvertebrate Communities to Hydrological and Water Quality Variations in a Heavily Regulated River

Aquatic macroinvertebrates are frequently used to evaluate river system conditions and restoration project performance. A better understanding of macroinvertebrate community responses to multiple stressors is a primary challenge for river science. In this paper, macroinvertebrate responses to hydrological and water quality variability were studied in the regulated Oglio River (northern Italy). We hypothesized that in regulated rivers the hydrological, rather than the physico-chemical conditions, would affect macroinvertebrate communities and biomonitoring tools (taxonomic metrics and functional indices). Repeated sampling (six times a year) was performed at four sites downstream of four dams in a 30 km river stretch during 2014 and 2015. Data were analysed using a linear mixed effect framework, to take into account random variation due to site and sampling date, and with multivariate analysis to track changes in community structure. A total of 69 families and 134,693 organisms were identified. The investigated metrics were mainly affected by the coefficient of variation of discharge, minimum discharge, ammonium, and temperature. The short-term dynamics of hydrological and physico-chemical variables were generally less important than the overall random effects as drivers of macroinvertebrate-based metrics. However, the relevance of a random effect (site, time, their interaction) differed depending on the biological metrics analysed. Understanding potential differences in response to short term and short stretch conditions would benefit biomonitoring and restoration procedures in both regulated and natural rivers.


Introduction
The hydrological alteration of river systems through impoundments and diversions represent one of the most important impacts on fluvial ecosystems [1] and a better understanding of the ecological responses is a primary challenge for river science [2]. More than 50% of large rivers worldwide are flow-regulated or dammed [1], and thousands of major dams are planned or under construction [3]. Projects on environmental flows, riverine restoration and those focussed on improving hydrographic connection, habitat heterogeneity and ecological integrity often centre on river hydrology and flow management [4][5][6]. However, possible relationships among hydrological, physico-chemical variables, and the community responses of river biota are yet to be fully understood [5,7,8]. Aquatic macroinvertebrates are model organisms for evaluating the ecological status of river ecosystems, as well as for detecting anthropogenic impacts of very different origin, spanning organic matter and nutrient contamination, deoxygenation and alteration of the flow regime [9][10][11][12][13]. Wide variations in river discharge and water velocity have been shown to alter the physical structure of the riverbed through either sediment clogging or flushing and to affect macroinvertebrate community structure [14,15]. Wide and sudden changes in river flow can dislodge benthic organisms thus triggering their drift and resulting in a local reduction of species richness and abundance [16]. Hydraulic force can also cause changes in food availability for benthic fauna which indirectly shapes community structure (e.g., abundance of benthic algae and fine particulate organic matter selecting different trophic functional groups) [17]. Together with the alteration of hydraulic features of watercourses, flow alteration may also affect the physico-chemical characteristics of water. For instance, changes in oxygen content and thermal regime have been reported by many authors [18][19][20] with the magnitude of change depending on basin size and dam typology (e.g., water released from the bottom or top of the dam).
In this work, we aimed to assess the effect of hydrology and water quality on macroinvertebrate community and biomonitoring tools along a regulated river system. In order to have a complete overview of the macroinvertebrate community, the taxonomic and functional community characteristics were considered. Numerous taxonomic metrics widely used in biomonitoring programmes and functional indices (functional richness and redundancy) were explored allowing us to go beyond a traditional taxonomic approach. Functional measures at the community level are related to ecosystem functioning and stability [21] and can reveal community responses to a variety of stressors [22,23]. Specifically, we aimed to study and test the different relevance of multiple stressors (hydrological, chemical, physical variability) on (i) community metrics and indices and (ii) community structure. We hypothesise that the flow/hydrological variables were the major drivers affecting metrics used to describe the macroinvertebrate community compared to physico-chemical characteristics of water. We expected the functional richness response to be concordant with the other taxonomic metrics (at least in relation to taxonomic richness following [24]). In addition, we postulated that functional redundancy can provide complementary information, and anticipate that it will increase with discharge, likely due to habitat homogenization. Moreover, as the sites are subject to similar hydro-morphological alterations and are located along the same river stretch within the same hydro-ecoregion for biomonitoring purpose, they would be expected to have similar macroinvertebrate communities. Indeed, disentangling the effect of multiple stressors on a macroinvertebrate community is crucial to provide water managers and river basin authorities with reliable information for conserving the integrity of watercourses and to plan long term ecological restoration. Understanding potential differences in responses to short-term and short-stretch conditions would also benefit biomonitoring approaches in both regulated and natural rivers.

Study Area
The study was performed along the Oglio River (Lombardy region, Northern Italy), a 280 km long tributary of the Po River. Along its course Lake Iseo originates, the sixth Italian lake in size (65.3 km 2 , 25 km long). This study was carried out downstream of the lake outlet that hosts six hydroelectric generation plants. The average annual production of these plants sums to 14.5 GWh. The oldest hydroelectric generation plant was built in 1933, while the most recent dates to 1984. More details Water 2019, 11, 1478 3 of 18 on the entire Oglio river basin and on the stretches downstream to the lake outlet are available in Guareschi et al. 2014 [24].
The study was performed at four sites each located downstream of different hydroelectric generation plants (Figure 1). These sites were located in the hydro-ecoregion called "Po Plain" following the Italian government designation (Italian Ministry Decree 260/2010), at an altitude ranging from 181 m (site 1) to 128 m (site 4). The wet riverbed in each site resembles that of the active river channel at the same site with a wide range of discharges due to morphological constraints, e.g., riverbed deepening and the river flowing through a single channel. The width of the riverbed was nearly 30 m at all investigated. The riverbed grain/clast size has been reported according to the Italian standard for biomonitoring, using the intermediate axis length-L [25]. Three sites were composed of coarse-grained materials classified as megalithal (L > 40 cm), macrolithal (20 cm < L < 40 cm) and mesolithal (6 cm < L < 20 cm), while the substrate of the fourth site was mainly mesolithal, microlithal (2 cm < L < 6 cm) and cobble (0.2 cm < L < 2 cm). The percentage of other mineral (e.g., sand, silt) or biotic (e.g., macrophytes, CPOM) substrates were less than 5%, thus limiting the substrate patchiness and its potential effect on macroinvertebrate community composition. Riparian vegetation was present at both side of the riverbed at each site, with the exception of site 1 that had riparian vegetation only on the left side. Considering the width of the active channel, the riparian vegetation cover was less than 10% at each site. The stream slope in the studied reaches was in all case <2.5% and all the sites were located within the first 1.5 km downstream of infrastructure. Agriculture is the main land use in the area and agricultural land represents about 60% of the basin area downstream of Iseo Lake [26] with higher values at the last three studied sites.  [24]. The study was performed at four sites each located downstream of different hydroelectric generation plants (Figure 1). These sites were located in the hydro-ecoregion called "Po Plain" following the Italian government designation (Italian Ministry Decree 260/2010), at an altitude ranging from 181 m (site 1) to 128 m (site 4). The wet riverbed in each site resembles that of the active river channel at the same site with a wide range of discharges due to morphological constraints, e.g., riverbed deepening and the river flowing through a single channel. The width of the riverbed was nearly 30 m at all investigated. The riverbed grain/clast size has been reported according to the Italian standard for biomonitoring, using the intermediate axis length-L [25]. Three sites were composed of coarse-grained materials classified as megalithal (L > 40 cm), macrolithal (20 cm < L < 40 cm) and mesolithal (6 cm < L < 20 cm), while the substrate of the fourth site was mainly mesolithal, microlithal (2 cm < L < 6 cm) and cobble (0.2 cm < L < 2 cm). The percentage of other mineral (e.g., sand, silt) or biotic (e.g., macrophytes, CPOM) substrates were less than 5%, thus limiting the substrate patchiness and its potential effect on macroinvertebrate community composition. Riparian vegetation was present at both side of the riverbed at each site, with the exception of site 1 that had riparian vegetation only on the left side. Considering the width of the active channel, the riparian vegetation cover was less than 10% at each site. The stream slope in the studied reaches was in all case <2.5% and all the sites were located within the first 1.5 km downstream of infrastructure. Agriculture is the main land use in the area and agricultural land represents about 60% of the basin area downstream of Iseo Lake [26] with higher values at the last three studied sites.

Sampling Activity: Biological, Hydrological and Water Quality Data
To fulfil the objectives of the research we planned an intensive biological sampling design that involved collecting 480 macroinvertebrate samples in 12 sampling campaigns. Macroinvertebrates were collected at investigated sites on an approximately monthly basis, on six dates from April to September during 2014 and 2015 (n = 12). A Surber net with a 0.05-m 2 area was used to collect macroinvertebrate samples. Ten replicates were randomly collected per site during each sampling date at least 2 m apart to minimise spatial autocorrelation [27]. Each replicate was stored separate from one another, and macroinvertebrate samples and associated material were preserved in 70% ethanol to be examined under a stereoscope in the laboratory. All the macroinvertebrate individuals were identified to the family level, except for Hydracarina, following the taxonomic classification proposed by Tachet [28]. This taxonomic level is the official prerequisite for riverine biomonitoring purposes in Italy and is commonly used in Europe and elsewhere (see [25,29]). Moreover, a good

Sampling Activity: Biological, Hydrological and Water Quality Data
To fulfil the objectives of the research we planned an intensive biological sampling design that involved collecting 480 macroinvertebrate samples in 12 sampling campaigns. Macroinvertebrates were collected at investigated sites on an approximately monthly basis, on six dates from April to September during 2014 and 2015 (n = 12). A Surber net with a 0.05-m 2 area was used to collect macroinvertebrate samples. Ten replicates were randomly collected per site during each sampling date at least 2 m apart to minimise spatial autocorrelation [27]. Each replicate was stored separate from one another, and macroinvertebrate samples and associated material were preserved in 70% ethanol to be examined under a stereoscope in the laboratory. All the macroinvertebrate individuals were identified to the family level, except for Hydracarina, following the taxonomic classification proposed by Tachet [28]. This taxonomic level is the official prerequisite for riverine biomonitoring purposes in Italy and is commonly used in Europe and elsewhere (see [25,29]). Moreover, a good relationship of higher taxonomic levels (e.g., family) with environmental variables has been observed in ecohydrological and ecological studies of river ecosystems [11,30,31].
All data related to discharges were provided by the Consorzio dell'Oglio, the official Basin authority that regulates discharges from Lake Iseo. Data and details on water regulation can be observed in Figure S1, where the natural and regulated water levels at the Sarnico dam (Lake Iseo outlet), together with the water inflow and outflow from the lake, are reported for hydrological years 2013-2014 and 2014-2015. The discharge descriptors used were (i) mean (Q mean ) (ii) maximum (Q max ) (iii) minimum (Q min ) and (iv) coefficient of variation (Q cv ).
Dissolved oxygen (% O 2 ), pH, temperature, and electrical conductivity were determined in situ using a multiparameter probe (YSI 556 MPS, Yellow Springs, OH, USA). Four replicate water samples were collected concurrently with macroinvertebrate sample collection. Water samples were filtered in situ (glass fiber filters GF/F Whatman, Maidstone, England, porosity 0.7 µm), stored refrigerated and processed within 6-8 h after sampling. Ammonium (NH 4 + ) and nitrate (NO 3 − ) were subsequently determined on filtered water while total phosphorus (TP) was determined on unfiltered water samples [32,33].

Macroinvertebrate Metrics and Indices: Taxonomic and Functional Approaches
Total family richness, abundance, ASPT (Average Score Per Taxon) and EPT family richness (Ephemeroptera, Plecoptera, Trichoptera family richness) were calculated. All these metrics represent single tools that frequently compose more complex multi-metric indices used in freshwater biomonitoring to fulfil the European Water Framework Directive 2000/60/EC (e.g., [34][35][36]). The ASPT represents a sensitive taxa index that is calculated for each sample based on the BMWP (Biological Monitoring Working Party) divided by the number of scoring families [37], it is a core metric for the Italian multimetric STAR_ICMi index [36] but also represent an independent index in Spain [38] and UK [39] and it is frequently used to indicate organic pollution.
Functional responses were also investigated through functional richness (FRich) and functional redundancy (FRed) indices. FRich is defined as the amount of functional space filled by the community [40], while FRed as the difference between taxonomic diversity and functional diversity [41] and it is related to ecosystem stability, resistance and resilience [42]. FRich and FRed were characterized using traits proposed by Tachet et al. [28]. However, to calculate FRed we selected five "effect traits" (size, dispersal, locomotion, food and feeding habits) based on Schmera et al. and Hevia et al. [43,44]. Functional indices calculation was based on the community-level weighted means of trait values (CWM). CWM matrix with the proportion of each trait character in each sampling site was obtained by crossing the "Taxon × traits" and "taxon × site" matrices. The relevance of the functional approach has been already stressed in some biomonitoring indices (e.g., France, [35]) in which macroinvertebrate trait metrics are directly integrated. Moreover, functional redundancy has recently been proposed as a biomonitoring tool to monitor responses in the riparian vegetation [22] and tested on macroinvertebrate communities in a highly regulated Mediterranean basin [45].

Statistical Analysis
All the taxonomic metrics and functional indices considered were analysed as a function of discharge and physico-chemical predictors within a linear mixed model framework to avoid pseudoreplication issues (LMM; [46]). By taking into account the characteristics of the experimental design, sampling date, site and their interaction were considered random effects, while predictors were considered as fixed. Prior to the analyses, collinear covariates were removed by stepwise selection based on the variance inflation factor (VIF). Abundance was log-transformed and fixed effects were scaled prior to analysis. Bayesian inference based on Markov chain Monte Carlo (MCMC) was used to carry out the analysis because it (i) allows an easier handling of random effects components, and (ii) provides full uncertainty propagation. Normal priors where used for the fixed effect and the inverse Wishart for random effect. A single chain was used and 1 × 10 8 iterations were carried out, with a burn in of 5000 and a thinning interval of 500. Such and high number of iterations was used to assure convergence of random effects. The convergence of the chain was assessed by visually inspecting the traceplot, using the Heidelberger and Welch's convergence diagnostic and by examining to the effective sample size estimates. Results are expressed as 95% credible intervals that represent the interval where there is a 95% chance that the true value of the parameter will be within it [47]. In our work credible intervals that not overlap 0 were considered significant.
Conditional and marginal coefficients of determination (pseudo R 2 ) were calculated according to Nakagawa and Schielzeth [48] with the aim to partition the variance explained by environmental (fixed) and latent (random) variables. Briefly, marginal R 2 represents the variance explained by fixed effects, while conditional R 2 represents the variance explained by both fixed and random effects. The variance explained by random effects was further partitioned into components related to site, time and the time/site interaction.
A PERMANOVA using Bray-Curtis dissimilarities was performed to test if the site, time, and their interaction significantly affect macroinvertebrate community structure. Non-metric multidimensional scaling (NMDS) was performed to identify patterns in community composition. Vectors of the covariates were fitted to the NMDS and performed separately for each site by using the function envfit in the vegan package with the aim of assessing the importance of each covariate in structuring macroinvertebrate community at the local scale. The contribution of each covariate was assessed using the square correlation coefficient (r 2 ).
Finally, an indicator value analysis (IndVal) was carried out to identify indicator families of the investigated sites [49]. The indicator value for abundance data is the product between two quantities, where the first is the mean abundance of the families at the target site group divided by the sum of the mean abundance values from all the groups, and the second is the relative frequency of occurrence of the family at the target site [50]. The IndVal analysis was run with the data aggregated per time and site. However, the presence of an indicator family at a site does not necessarily mean that this family is absent from the other sites, but that its abundance and detection frequency at that identified site are higher.
All the statistical analyses were performed using the statistical computing software R [51]. VIF selection were performed with the package "usdm" [52]. Bayesian LMMs were performed with the package "MCMCglmm" [53] and the package "coda" was used to assess models convergence [54]. Package "vegan" [55] was used for PERMANOVA, nMDS and vector fitting while "indicspecies" [50] was used for the IndVal analysis. Package "ggplot2" was used to plot the results [56]. R scripts used to calculate EPT, ASPT, FRic and FRed are available at https://github.com/alexology/biomonitoR.  Table 1. Water was well oxygenated with an alkaline pH and low conductivity (  Q max , Q mean and electrical conductivity were excluded from the statistical analysis following the results of the stepwise selection based on VIF. After excluding these descriptors VIF results were lower than 2 and the maximum and minimum correlations among predictors were −0.43 (T-NO 3 − ) and 0.003 (pH-NO 3 − ).  Table S1. Among the EPT taxa, Plecoptera were only represented by the family Leuctridae, which consists of species characterised by very different ecological requirements mostly not particularly pollution-sensitive [57]. Heptageniidae, Leptophlebiidae and Potamanthidae were the Ephemeroptera families that exhibited high sensitivity to anthropogenic disturbance based on BMWP scores. However, the last two families were present at very low abundances (<15 individuals) and only on a few sampling occasions.

Macroinvertebrate Community: Responses, Ordinations and Indicator Taxa
LMM results are reported in Figure 2. Credible intervals of hydrological covariates were differed from 0 for all the metrics and indices tested (Q cv and Q min ). Among water quality variables, credible intervals differed from 0 for NH 4 + (all the metrics and indices except FRich and FRed) and T (only for abundance and FRed). Metrics and indices were in general negatively related to covariates, despite this was not true for FRed. For example, a negative relationship with the hydrological variable Q min means that an increase in minimum discharge causes a decrease of the target metric or index. The R 2 values, reported in Table 3, show the importance of the random effects (site, time and the time/site interaction) compared to the fixed effects (discharge, physico-chemical predictors) to explain the variation among the dependent variables. Moreover, the relevance of a random effect depends on the biological metrics analysed (Table 3).  Site, time and their interaction explained 9%, 23% and 18% of the total variance in the PERMANOVA analysis, respectively. Community structure at the investigated sites differs in their temporal trajectories despite having similar environmental characteristics. NMDS was fitted on 3 axes with a stress value of 19%. The relevance of site and time on community structures can be also observed in Figure 3 with a remarkable vertical shift along the time gradient being observed from  Site, time and their interaction explained 9%, 23% and 18% of the total variance in the PERMANOVA analysis, respectively. Community structure at the investigated sites differs in their temporal trajectories despite having similar environmental characteristics. NMDS was fitted on 3 axes with a stress value of 19%. The relevance of site and time on community structures can be also observed in Figure 3 with a remarkable vertical shift along the time gradient being observed from time 1 to 6 and from time 7 to 12 (Figure 3b,c). Vector fitting performed separately for each site highlighted the importance of site-specific characteristics in structuring the macroinvertebrate community (Table 4). NO 3 − , NH 4 + and temperature strongly affected macroinvertebrate community structure at all the investigated sites while hydrology was important for Site 1 (Table 4), the closest to the outlet of Lake Iseo. The IndVal analysis selected a total of 23 taxa with significant indicator value (IV) > 0.5, from just one (at site 2) to ten (at site 4) (complete details in Table 5). The indicator family group highlighted for site 1 was characterised by the presence of lentic families (e.g., Dreissenidae, Planorbidae, Acroloxidae, Hydridae), while site four taxa were largely associated with lotic (e.g., Leuctridae, Heptageniidae) ( Table 5). At site 2, one indicator taxon was identified (Psychomyiidae) while site 3 presented a heterogeneous list of five indicator families, in both cases with lower IV in most of the cases.
highlighted the importance of site-specific characteristics in structuring the macroinvertebrate community (Table 4). NO3 − , NH4 + and temperature strongly affected macroinvertebrate community structure at all the investigated sites while hydrology was important for Site 1 (Table 4), the closest to the outlet of Lake Iseo.  2015 (c). A confidence level of 0.5 was adopted for the ellipsoid to increase the readability of the plot and to better assess spatial and temporal trends.
The IndVal analysis selected a total of 23 taxa with significant indicator value (IV) > 0.5, from just one (at site 2) to ten (at site 4) (complete details in Table 5). The indicator family group highlighted for site 1 was characterised by the presence of lentic families (e.g., Dreissenidae, Planorbidae, Acroloxidae, Hydridae), while site four taxa were largely associated with lotic (e.g., Leuctridae, Heptageniidae) ( Table 5). At site 2, one indicator taxon was identified (Psychomyiidae) while site 3  2015 (c). A confidence level of 0.5 was adopted for the ellipsoid to increase the readability of the plot and to better assess spatial and temporal trends.

Macroinvertebrate Responses: The Importance of Fixed and Random Effects
Many factors act on different spatial and temporal scales to shape macroinvertebrate communities in running waters. Within this framework, our research has demonstrated that the short-term dynamics of hydrological and physico-chemical variables seem to be rarely the most important drivers of the macroinvertebrate community. In contrast, we found that site, time, and site/time interactions (the so-called random effects) explained a relevant proportion of variance especially for biomonitoring tools. However, the relevance of a random effect (site, time, their interaction) differed depending on the biological metrics analysed (Table 3). For example, ASPT and EPT richness, the metrics most influenced by random effects, are predominately affected by the "site" (and poorly related to "time") while its effect seems extremely weak for abundance values and FRich (Table 3). This may suggest that metrics characterise spatial or temporal pattern of macroinvertebrate community differently and they need to be carefully chosen based on biomonitoring procedures and restoration targets (see below).
However, based on all the predictors tested (fixed effects), the role of selected hydrological variables on biomonitoring tools can be highlighted. We found partial agreement with our initial hypothesis, the coefficient of variation (Q cv ) and the minimum values of discharge (Q min ) were negatively related to all metrics and indices considered, except for FRed that displayed a contrasting pattern.
Responses of macroinvertebrate communities to hydrological factors is reported to be highly variable, spanning either increases or decreases in both abundance and diversity [58]. In our research on regulated river reaches, the negative effect of hydrological variables could be partially ascribed to periods of high flow that occurred mainly in the very wet year 2014. The detrimental effect of floods was better represented by Q cv and elevated Q min values, according to Konrad et al. [59]. The increase in shear stress during a flood is expected to cause macroinvertebrates to drift from the benthic zone into to the water column [60]. Dislodged organisms may be transported downstream, with a reduction in their abundance at a given site, and with potential effects on taxa richness in some instances (e.g., [61]). The sudden increase in near-bed velocity increases both the probability of dislodging some lentic macroinvertebrates from the substrate [62] and the abundance of filter-feeders adapted to high water velocity [63]. Changes in discharge (probably related with high values of Q cv ) may modify the availability of oviposition sites for aquatic insects (e.g., emergent rocks) altering egg supply and subsequently the biological communities [64]. Q cv could be related to improvements in habitat heterogeneity and formation (due to dynamic and variable discharges), although it could also reflect the dramatic changes in discharges that negatively affect the invertebrate communities in heavily regulated rivers. In this context, defining and characterizing this 'natural dynamic flow regime', which is crucial for maintaining biodiversity and ecological processes, into quantitative flow descriptors remains a challenge for river science (see [5] and reference inside).
The reduced abundance following the flood was not surprising, and has been reported by other authors [65,66]. In our study area sensitive taxa (EPT richness and ASPT) decreased with increasing Q cv and Q min . However, Suren and Jowett [67] found EPT metrics to be resilient to flood alterations in a New Zealand river. The decrease we observed could be partially affected by our experimental design, in which we considered the mean response of the dependent variable related to an area of 0.05 m 2 . A rarefaction of sensitive taxa at 0.05 m 2 might not match a decrease in sensitive taxa at a larger scale. However, similar patterns between both metrics were expected as EPT taxa can affect ASPT scores as most have higher BMWP scores. As expected, the response of FRich was largely concordant with taxonomic metrics, while in general, a contrasting pattern was observed for FRed. Community functional redundancy quantifies the degree of niche overlap in functional space and reductions have been highlighted along stress gradients (e.g., [68]). In the Oglio River, increased values of Q cv and Q min appeared to produce a net increase in niche overlap among taxa. In response to our hypothesis, increasing levels of Q min in regulated rivers could be related to habitat homogenization or stabilization with few possibilities for new and complementary taxa with different functional features to enter sites with consequent effects on FRed. In contrast to FRich, Q cv had a positive effect on FRed. This may be related to abrupt changes that select and filter specific and resistant taxa progressively. Concordance between FRich and FRed in aquatic communities have been observed when responding to different stresses [68] while in our research these functional indices behaved differently depending on the variables tested.
In agreement with our findings, Bruno et al. [69] also found different responses of FRich and FRed when dealing with flow regulation in riparian vegetation communities. Furthermore, in our study FRed was affected more by random effects than functional richness, the "site" being particularly important. The negative relationship between temperature and aquatic faunal abundance is probably due to the emergence of insects from April to September in both years. The positive relationship between temperature and FRed could be explained by the distribution of common and stress-resistant taxa. Possible effects of seasonal habitat variability on both taxonomic and biological trait composition have been already stressed by other authors in other systems, especially for abundance [70].
For taxonomic metrics, hydrological variables and NH 4 + provided similar significant patterns. NH 4 + significantly reduced abundance, family richness, EPT richness and ASPT while no relationships were observed for functional indices. A reduction in macroinvertebrate richness and sensitive taxa could be attributed to the toxicity exerted by NH 4 + on aquatic organisms [71,72]. However, this is unlike in the Oglio River where NH 4 + concentrations were much lower than those reported in literature [73,74].
The river reach studied in this research was fed by the surface waters of Lake Iseo, in which NH 4 + is often undetectable [75]. Despite a toxic effect not being excluded, the relationship between NH 4 + and the macroinvertebrate metrics studied probably reflects a latent ongoing process that is strictly related to NH 4 + , as suggested by Friberg et al. [10] and Beketov and Liess [76] for other rivers. Such ongoing processes are probably related to the influence of minor polluted tributaries, draining heavily exploited farmland [26], and partially untreated sewage from overflows during heavy rain events [77].

Differences in the Similarity: Cross-Site Comparison
Our results highlighted that even in short river stretches (~30 km) of heavily regulated rivers, the benthic macroinvertebrate community responded differently depending on local characteristics, time, and the interaction of these factors. The importance of covariates depends on site-specific characteristics as demonstrated by the envfit analysis. Numerous findings emphasized these differences among sites. Community composition showed spatial and temporal trends, mainly related to NO 3 concentration and temperature, which can be proxies of both river-groundwater interactions and the time of the year (e.g., the lower the base flow, the higher the input of groundwater). As for NH 4 + , a direct effect of NO 3 − on macroinvertebrates should be excluded given the low concentrations recorded at the investigated sites. However, the influence of NO 3 − on macroinvertebrate community composition is especially evident at sites further away from Lake Iseo (e.g., site 3, Table 4) as previously demonstrated by Laini et al. [36]. Groundwater inputs have been reported to support macrophyte assemblages in the same system by providing habitat heterogeneity [8,78], these may in turn indirectly affect macroinvertebrate community structure [79]. According to Bonada et al. [80], following Dufrêne and Legendre [49], IV ≥ 0.5 (in case of "indicespecies" package approach) can be considered as a threshold to judge adequate an indicator taxon. Most of them come from site 1 (seven taxa) and site 4 (10 taxa) testifying their differences in community composition. These results, along with PERMANOVA outputs (significant difference among sites communities) and NMDS ordinations (Figure 3, Tables 4 and 5), support the conclusion that there are differences among sites subject to similar pressures in the same river stretch. The drift of the organisms from Lake Iseo affected the community structure at site 1, as clearly stressed by the IndVal analysis, due to the presence of taxa characteristic of lentic water, such as Dreissenidae, Planorbidae and Acroloxidae (and even Hydra). This pattern can be interpreted as a mass effect [81], and similar findings for Dreissena polymorpha have been recorded downstream dams in some Croatian rivers related to the reproduction stage of planktonic larvae in reservoirs [82].
The presence of Lake Iseo, moreover, prevents the drift of lotic taxa from the upstream rivers to the investigated sites and indirectly provides the river with lentic taxa. Numerous lotic taxa appeared at just site 4 alone, the farthest from the lake, similar to the pattern reported by Holt et al. [83] in the Chattahoochee River (Georgia, USA) where no EPT indicator taxa existed directly below the dam but increased with the distance downstream.
These results, together with the importance of the site/time interaction in LMM, generally suggested that local dynamics play a primary role in shaping macroinvertebrate communities, even at sites with similar hydrological and physico-chemical features as well as near-uniform substrate conditions, climate and altitude. Macroinvertebrate taxa richness and density have been recognised to be highly influenced by local and small scale spatial features in different systems and geographical areas [84][85][86][87].

Remarks and Implications for River Monitoring and Management
To contrast anthropogenic pressures on freshwater ecosystems, biomonitoring and ecological restoration are essential tasks for sustainable water resource management, aquatic biodiversity conservation and environmental legislation. Multiple natural and anthropogenic factors shape the structure of macroinvertebrate community in the flow regulated Oglio River. Despite the contribution of hydrological and physico-chemical variables cannot be easily disentangled the outcomes of this study stress the significance of hydrological variables for all biomonitoring metrics, both taxonomic and functional.
Among the studied metrics, Family richness, EPT richness and ASPT displayed complete concordance to all stressors tested in the study, illustrating how hydrological variables and NH 4 + as important predictors of macroinvertebrate community. The response of FRed stands out as it displayed the opposite trend compared to all other metrics and indices analysed. This lack of concordance indicates that caution is required when using it as a surrogate for other parameters, when assessing river quality or restoration project outcomes. It must be acknowledged that our research was performed at a high taxonomic resolution (family level). When available, species information should be used, however our findings at family level for both taxonomic and functional measures indicate common patterns and reliable results have been achieved, at least in the context of rapid bioassessment.
Improved understanding of spatial variability is a scientific challenge in environmental flows and water resource management [88]. Even in a short river section with near-uniform substrate conditions, climate and altitude, the macroinvertebrate community responds differently to environmental factors depending on very local river characteristics. So called random effects (site and time, here) exerted strong effects despite with some difference depending on the metrics considered. The results indicate that in biomonitoring procedures and environmental impact assessment, local conditions (e.g., small scale variability) should be considered given the difficulty of generalizing univocal responses, at least, along the Oglio River. Further research is required to disentangle the effect of multiple stressors on macroinvertebrate community and to provide reliable information and tools for preserving the integrity of running water ecosystems. Therefore, local characteristics and small-scale variables should be integrated in basin scale approaches when undertaking restoration projects, flow managements and when setting the bases for the longitudinal delimitation of heavily modified water bodies (HMWB, sensu WFD 2000/60/CE) downstream of reservoirs.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/7/1478/s1, Figure S1: Water influx to and water outflow from the lake Iseo and regulated and natural water levels at Sarnico dam, Table S1: Taxa collected at the investigated sites.