Climate Change and Water Exploitation as Co-Impact Sources on River Benthic Macroinvertebrates

: Climate change can affect freshwater communities superimposing on other major stressors, such as water exploitation, with effects still poorly understood. The exacerbation of naturally-occurring periods of low ﬂows has been reported as a major hydrological effect of water diversions, with severe impacts on river benthic macroinvertebrate communities. This study aimed at assessing long-term modiﬁcations of low-ﬂow events in a large lowland Italian river possibly associated to climate change and the effects of these events, intensiﬁed by water withdrawals, on benthic macroinvertebrates. A 77-year dataset on daily discharge was thus analyzed through Mann-Kendall test and Sen’s method to investigate modiﬁcations of the main hydrological parameters. Moreover, macroinvertebrates were collected during the low-ﬂow periods that occurred from 2010 to 2015 at three sites downstream of water withdrawals, representing three different conditions of hydrological impairment. After assessing possible differences in taxonomical and functional composition between sites and impairment conditions, redundancy analysis and ordinary least squares regression were performed to link benthos metrics to environmental (hydrological and physico-chemical) characteristics. An increase in the duration of the low-ﬂow periods and reduced summer ﬂows were detected on the long term, and the magnitude of low ﬂows was signiﬁcantly altered by water withdrawals. These hydrological features shaped both structural and functional characteristics of benthic assemblages, highlighting the need for a more environmentally-sustainable water resource management in the current context of climate change.


Introduction
Biological communities of lotic ecosystems are subjected to numerous anthropogenic pressures, including pollution, introduction of invasive alien species, morphological and hydrological alterations, and climate change. Particularly, the last two are strictly interconnected, since modification of temperatures, precipitations, and flow regimes induced by climate change can affect both water availability and freshwater demand for human activities [1]. Low flows are one main hydrological components that can be significantly influenced by both water exploitation and climate change. Low-flow periods, which naturally occur in flow regimes, are magnified by off-stream water diversion, particularly in irrigated areas where low flows naturally take place in summer, as it commonly occurs in the Mediterranean region [2]. Despite modelling uncertainty in hydro-climatic scenarios, there is general consensus about expected exacerbation of low-flow periods over the next decades in central and southern Europe, e.g., [3][4][5][6][7][8].
From an ecological point of view, naturally occurring low flows act as an environmental filter by selecting species owing traits that make them able to survive to such extreme

Study Area
The Ticino River is the main tributary of the Po River (northern Italy) in terms of flow rate (MANF at Po confluence: 348 m 3 s −1 ). It flows from the Swiss Alps to Lake Maggiore for~90 km and then from the lake to the Po River for~110 km ( Figure 1). Its watershed area is 8172 km 2 , including mountainous areas from the source to the lake outlet and lowland agricultural areas. Urbanized areas are mainly located near the lake and on the east side of the lowland course. Climate in the Ticino River watershed ranges from cold (Df, in the upland catchment) to temperate (Cf, in the lowland catchment, including the study reach), based on Köppen climate classification [26].
Since 1943, the lake outflow is regulated by the Miorina Dam ( Figure 1) to optimize the water use in the downstream section of the river; no water diversion takes place at this dam. Water releases through the dam depend on the water level of the lake, according to limits established by laws. Discharge is controlled through movable gates that are fully opened during floods. The mean annual lake outflow is 280 m 3 s −1 .
relevant diversion (maximum 44 m 3 s -1 for mixed irrigation and hydropower use) occurs 22 km downstream of the Panperduto Dam, by the Langosco Canal intake (active since 1665; Figure 1). Mandatory monthly-modulated MFs are currently released below both the Panperduto Dam and the Langosco Canal intake. However, streamflow generally exceeds MFs during the periods of high water availability, i.e., mainly in spring due to snowmelt and in autumn due to rainfall. Conversely, two periods of low flows usually occur in winter and in summer. In the lowland course (i.e., downstream of Lake Maggiore), the authorized maximum cumulative water diversion approximately equals the mean annual outflow of the lake. Specifically, about 4 km below the Miorina Dam water is diverted to the western side of Ticino River (max. 70 m 3 s −1 ) through the intake of the Regina Elena Canal, mostly for irrigated agriculture. About 3 km downstream, the largest off-stream diversion takes place at the Panperduto Dam (Figure 1), which has been active since 1884, and currently allows for a maximum diversion to the eastern side of Ticino River of approximately 60% of the river mean annual discharge (i.e., 170 m 3 s −1 ), for agriculture and hydropower. A further relevant diversion (maximum 44 m 3 s −1 for mixed irrigation and hydropower use) occurs 22 km downstream of the Panperduto Dam, by the Langosco Canal intake (active since 1665; Figure 1). Mandatory monthly-modulated MFs are currently released below both the Panperduto Dam and the Langosco Canal intake. However, streamflow generally exceeds MFs during the periods of high water availability, i.e., mainly in spring due to snowmelt and in autumn due to rainfall. Conversely, two periods of low flows usually occur in winter and in summer.
We selected three sampling sites ( Figure 1)  amounting to 2-4% of the MANF is released, which is gradually increased through further downstream contributions, reaching the same MF values as in MF1 in a few km.
River morphology was instead comparable among sites. Except for the presence of minor riprap areas, the riverbed morphology is predominantly natural in all sites, with alternation of riffles, runs, and pools, and the substrate is mainly composed of cobbles (6-20 cm diameter).

Data Collection
At the three study sites, benthic macroinvertebrate sampling was carried out during low-flow periods, on one to three occasions according to the length of the period, from 2010 to 2015 (N = 20 samples per site, see Table S1). Benthic macroinvertebrates were collected with a Surber sampler of 0.05 m 2 area and 500 µm mesh, following a quantitative multi-habitat protocol [27]. The adopted sampling method requires three steps: (1) visual quantification of the different substrates within a representative river reach, (2) proportional assignments of ten subsamples according to the areal distribution of the substrates, and (3) integration of the ten subsamples into a single sample, corresponding to a 0.5 m 2 total sampled area. The samples were preserved in formalin (4%) and analyzed in the laboratory. Individuals belonging to Plecoptera, Ephemeroptera, Turbellaria, and Hirudinea were identified to the genus level and the others to the family level using dichotomous keys [28] and photographic atlas [29].
On each benthos sampling date, temperature (T), dissolved oxygen (DO) concentration and saturation, pH, and electrical conductivity (EC) were measured in situ (at a single point per sampling site) using the YSI Professional Plus portable multiparameter probe (YSI, Yellow Springs, USA). Additionally, water samples were collected and analyzed in the lab according to standard methods [30] to measure the concentration of ammonia nitrogen (NH 4 + -N), nitrate nitrogen (NO 3 − −N), total nitrogen (TN), total phosphorus (TP), chemical oxygen demand (COD), and biochemical oxygen demand in five days (BOD 5 ).
Mean daily streamflow downstream of the Miorina Dam (i.e., the lake outflow, Figure 1) and at the three sampling sites were kindly provided by the Ticino Consortium for the periods 1943-2019 and 2010-2015, respectively. Specifically, the lake outflow and the discharge in the mentioned diversion canals were measured by stage gauging, starting in the 1940s. An additional gauging station located close to IF1 (Figure 1) was active during 2010-2015. Flows at the sampling sites MF1 and MF2 were then calculated by continuity.

Data Analysis
The streamflow time series (from 1943 to 2019) at the Miorina Dam was processed to detect temporal changes in the main hydrological parameters by using the software developed by The Nature Conservancy, which allows for the calculation of the indicators of hydrologic alteration (IHA, [31]). A preliminary check to detecting autocorrelation was performed by Ljung-Box test. Then, the Mann-Kendall test (or its modified version for significantly autocorrelated parameters [32]) and Sen's method were applied to test for significant temporal trends in the values of 47 hydrological metrics, i.e., the 33 IHA, the average monthly flows normalized by the average annual flow, and the frequency and duration of extreme low-flow periods (i.e., periods with flow values lower than the 10th percentile).
The IHA parameters were also calculated for the streamflows of the three sampling sites, for the six years of the macroinvertebrate study (from 2010 to 2015). A discriminant analysis (DA) of the calculated hydrological parameters was performed to assess possible differences between streamflow patterns in the three sampling sites. Moreover, possible differences of the water physico-chemical parameters measured concurrently to benthos sampling at the three study sites were evaluated using the Kruskal-Wallis test, followed by the Dunn test for pairwise comparisons. The structure of the zoobenthic communities was described through standard metrics, i.e., total density (N), family richness (S), and Shannon-Wiener diversity (H). We also considered the richness of families belonging to the orders Ephemeroptera, Plecoptera, and Trichoptera (EPT), the relative abundance of EPT taxa (%EPT), and the ratio of Chironomidae/Diptera individuals (C/D). The functional analysis was conducted by calculating the abundance of each functional feeding group (FFG) in the macroinvertebrate assemblages (active and passive filter feeders = AFF and PFF, gatherers-collectors = GAT, grazersscrapers = GRA, miners = MIN, predators = PRE, parasites = PAR, shredders = SHR, and xylophagous taxa = XYL), based on the information available in the freshwaterecology.info database v. 7.0-10/2016 [33]. Specifically, a fuzzy data coding was used, i.e., each family/genus in the database was associated to each of ten feeding types with a value ranging from 0 (no preference) to 10 (complete selectivity for one feeding type). Per each taxon, we multiplied this value by the number of individuals counted. The total number of individuals belonging to each feeding group was finally calculated summing up and dividing by 10 all the values referring to each feeding type. From the same database, information about voltinism was extracted to calculate the relative abundance of taxa with one or less generation cycle per year (i.e., semivoltines and univoltines) and that of taxa with more than one generation cycle per year (i.e., multivoltines). Moreover, based on the values of the above-mentioned functional groups, the following metrics proposed by Merritt et al. [34] were calculated: the ratio between primary consumers feeding on autochthonous primary producers and those feeding on allochthonous biomass, indicating the degree of autotrophy of the system (auto/heterotrophy-A/H); the ratio between secondary consumers and primary consumers, indicating control of predators on prey (top-down control-TD c ); and the ratio between the relative abundance of multivoltine organisms and that of semivoltines and univoltines, indicating the pioneering degree of the community (life cycle-LC).
In order to account for the variability due to sampling season and for seasonal differences in water uses, all samples were categorized into two sampling periods, as follows: samples collected from November to April, i.e., during cold months, when streamflow is mostly affected by withdrawals for hydropower (indicated as Non-Irrigation cold-NIcold-samples), and samples collected from May to October, i.e., during warm months, when large water withdrawals for irrigation take place (indicated as Irrigation warm-Iwarm-samples).
Compositional dissimilarity (both taxonomical and in FFGs) between macroinvertebrate samples was quantified using the Bray-Curtis index. This index ranges from 0 (indicating complete similarity) to 1 (indicating complete dissimilarity). Bray-Curtis distances were visually displayed in a two-dimensional ordination space using a non-metric multidimensional scaling (NMDS). Two-way analysis of similarities (ANOSIM) was performed to test the macroinvertebrate assemblages for significant differences in taxonomical and functional composition accounting for "site" and "period" factors. Similarity of percentages (SIMPER) was used to identify the taxa responsible for significant differences in the community structure.
Differences of the benthos metrics between sample groups (site*period) were analyzed by the Kruskal-Wallis test, followed by the Dunn test for pairwise comparisons.
For all the above-mentioned tests, significance level was set at p < 0.05 and Bonferroni correction of the p values was applied.
Environment-ecology relationships were investigated by redundancy analysis (RDA) and ordinary least squares (OLS) regression between benthos metrics and both hydrological and physico-chemical parameters. The latter comprised the above-mentioned spot measurements, carried out concurrently to benthos sampling. In the case of hydrological variables, twenty metrics (selected and modified from Schneider et al. [35]) were computed over the 90-days time-span before macroinvertebrate sampling, an adequate time interval for the detection of the effects of water diversion on macroinvertebrate assemblages [36]. Sorted by the five major categories of the hydrological regime (magnitude, rate of change, frequency, duration, and timing, [11]), the adopted flow metrics were: -Magnitude: mean flow (Q M ), coefficient of variation (Q CV ), minimum (Q min ), maximum (Q max ), 25th (Q P25 ), 50th (Q P50 ), and 75th (Q P75 ) flow percentiles, difference between maximum and minimum (∆Q), and mean flow of the sampling date (Q S ); -Rate of Change: mean and maximum increase (INC M and INC Max ) and decrease (DEC M and DEC Max ), and last increase (INC L ) and decrease (DEC L ) of flows between two consecutive days; -Frequency: number of low-flow days (FRE LF ), here defined as days with a flow lower than 10% of the MANF, and number of high-flow days (FRE HF ), defined as days with discharges at least three times larger than the median flow [37]; -Duration: maximum duration (number of days) of low flows (DUR LF-max ) and duration of the low-flow period immediately before the sampling (DUR LF-last ); -Timing: number of days from the last high-pulse (TIM HF ) day (this metric was not calculated over the previous 90 days but considering the whole dataset before each sampling date).
For the RDA, only two hydrological variables were selected, based on the related contribution to the two main axis of the principal component analysis (PCA-carried out with all the parameters) and the correlation with the other variables. Before RDA, all the variables were standardized.
For OLS regression, the best model approach according to Akaike information criterion was used. Unstandardized and untransformed data were used to estimate the coefficients directly in the original data units, and the model result was constrained to only one independent variable in order to find simple associations potentially informative for water resources managers and environmental protection authorities.

Long-term Analysis of Lake Maggiore Outflows
Considering the 77-year time series (1943-2019) of the Lake Maggiore outflows, seven out of the considered 47 hydrological metrics showed a significant (p < 0.05) monotonic trend (Table 1 and Figure 2). Table 1. Statistics of Ljung-Box, Mann-Kendall test (or its modified version [32]*), and Sen's method for the significant trends of the hydrological parameters concerning the Lake Maggiore outflows recorded during the period 1943-2019 (see  flows (by 26 days in 77 years). Accordingly, also the median duration of the extreme lowflow periods increased from 2 days in 1943 to 26 days in 2019. Moreover, a decrease of the mean monthly flow normalized by the annual mean was detected for September, when the monthly flow decreased on average from above to below the annual mean. The flow difference between consecutive days during flow decreases (fall rate) showed an increase from 12 to 21 m 3 s −1 , and the number of reversals decreased significantly from 98 in 1943 to 33 in 2019.

Figure 2.
Data and Sen's slope estimates of the hydrological parameters concerning the Lake Maggiore outflows recorded during the period 1943-2019 characterized by significant monotonic increasing or decreasing trend as detected by the Mann-Kendall test or its modified version [32] (see Table 1 Table 1 Only two metrics (i.e., Lo pulse # and reversals) were significantly autocorrelated according to Ljung-Box test ( Table 1).
The base flow index (i.e., the 7-day minimum flow divided by the mean annual flow) decreased from 0.42 to 0.31 according to the Sen's estimate. A significant decrease (from 6.28 to 0.98) in the frequency of low-flow periods (i.e., periods with discharge lower than the 25th percentile) was detected, along with an increase of the median duration of low flows (by 26 days in 77 years). Accordingly, also the median duration of the extreme low-flow periods increased from 2 days in 1943 to 26 days in 2019. Moreover, a decrease of the mean monthly flow normalized by the annual mean was detected for September, when the monthly flow decreased on average from above to below the annual mean. The flow difference between consecutive days during flow decreases (fall rate) showed an increase from 12 to 21 m 3 s −1 , and the number of reversals decreased significantly from 98 in 1943 to 33 in 2019.

Analysis of Low-Flow Periods Downstream of Water Withdrawals in the Ticino River
During the six investigated years (2010-2015), low flows in the Ticino River downstream of water withdrawals were mainly observed from January to April (winter-early spring) and from mid-July to October (summer-early autumn), consistently at the three study sites (Figure 3a). This pattern reflected that of natural flow variability but was sharpened in warm months due to maximum off-stream diversion, especially for irrigation ( Figure 3b).
The DA of the IHA confirmed the expected differences between the three sampling sites. Indicators related to the magnitude of low flows (particularly 3-day min and 1-day min) were the most important parameters differentiating IF1 from the other two sites. In fact, at IF1, these variables were higher than at MF1 and MF2 (i.e., on average, 18, 13, and 6 m 3 s −1 for 3-day min, and 17, 13, and 5 m 3 s −1 for 1-day min, respectively). In contrast, MF1 and MF2 were poorly discriminated, though at MF2 compared to MF1 (see Table S2  min, and rise and fall rates (on average, 10 vs. 25 m 3 s −1 for September mean flow, 7 vs. 25 m 3 s −1 , and 12 vs. 28 m 3 s −1 for rise and fall rates, respectively).
Differently from the IHA, most of the water physico-chemical parameters measured concurrently to benthos sampling during low-flow periods were not significantly different among the three study sites (Table S3). Only NO3 − −N was significantly lower at MF1 compared to the two downstream sites, and pH and % DO were significantly higher at MF2 than at IF1 (Kruskal-Wallis and Dunn test, p < 0.05). Differently from the IHA, most of the water physico-chemical parameters measured concurrently to benthos sampling during low-flow periods were not significantly different among the three study sites (Table S3). Only NO 3 − −N was significantly lower at MF1 compared to the two downstream sites, and pH and % DO were significantly higher at MF2 than at IF1 (Kruskal-Wallis and Dunn test, p < 0.05).

Analysis of Benthic Macroinvertebrate Communities during the Low-Flow Periods
NMDS did not reveal evident differences in taxonomical ( Figure 4a) and FFGs (Figure 4b) composition among the three sampling sites. However, benthic assemblages sampled during the irrigation period (i.e., from May to October, I-warm) showed a poorer variability compared to those sampled during the non-irrigation period (i.e., from November to April, NI-cold), both considering taxonomical and FFG composition; this occurred particularly at MF1, and, to a lesser extent, at MF2.
The two-way ANOSIM of taxonomical composition confirmed these patterns, indicating significant differences (R = 0.19, p = 0.0005) between sampling periods but not between sampling sites (R = 0.04, p = 0.07). This difference was mostly (>95%) related to the 20 most common families detected at the three sampling sites, as indicated by the SIMPER analysis (Table 2). Among these, some taxa (e.g., the Diptera Chironomidae and the Ephemeroptera Ephemerella) were more common in the NI-cold period, while others (e.g., the Ephemeroptera Baetis and the Plecoptera Leuctra) were more abundant during the I-warm period. cating significant differences (R = 0.19, p = 0.0005) between sampling periods but not between sampling sites (R = 0.04, p = 0.07). This difference was mostly (> 95%) related to the 20 most common families detected at the three sampling sites, as indicated by the SIMPER analysis (Table 2). Among these, some taxa (e.g., the Diptera Chironomidae and the Ephemeroptera Ephemerella) were more common in the NI-cold period, while others (e.g., the Ephemeroptera Baetis and the Plecoptera Leuctra) were more abundant during the Iwarm period.    Considering FFG composition, the two-way ANOSIM did not evidence significant difference both between sampling sites (R = 0.006, p = 0.37) and periods (R = 0.01, p = 0.32).
Among the biological metrics, the median values of S, H, and EPT decreased from IF1 (i.e., the river reach subjected to comparatively minor hydrological alteration) to MF2 (i.e., the river reach with the lowest MF values). Moreover, higher values of the mentioned metrics were observed in the I-warm period compared to the NI-cold one; however, these differences were not significant between sampling periods of fixed sites ( Figure 5). The PCA of the hydrological variables computed over the 90 days before benthos sampling and carried out separately for the two sampling periods, explained 65.7% of the overall variability in the datasets (see Table S4 and Figure S1). Since many variables were significantly correlated, only two parameters were selected for the RDA. Specifically, FRELF (i.e., the number of low-flow days) considerably contributed to PC1 and was significantly correlated to most of the other variables, and Qmin (i.e., the minimum flow value) considerably contributed to PC2. Both parameters characterized low flows at the study sites.
The RDA demonstrated that the variation of benthos metrics can be explained by physico-chemical parameters and by the considered hydrological variables in both sampling periods (I-warm = 63%, and NI-cold = 77%, Figure 6). A positive correlation between Qmin and H was detected as well as between nutrient load (mainly NO3 − −N) and LC and between the organic load expressed by BOD5 and TDc. Some differences between the two sampling periods could also be observed. For instance, in the I-warm period, %EPT was negatively correlated with NO3 − −N while in the NI-cold period it was positively correlated with BOD5 and Qmin.  % EPT in the I-warm period was higher than in the NI-cold period at MF1 and MF2 (with significant difference only for MF1, Kruskal-Wallis and Dunn test, p < 0.05), while a reversed pattern was detected at IF1. Conversely, LC was higher in the NI-cold than in the I-warm period at MF1 and MF2, and again the opposite was true for IF1, even if these differences were not significant ( Figure 5).

I-warm NI-cold I-warm NI-cold I-warm NI-cold
From the comparison of sites within the same sampling period, only H was significantly different (Kruskal-Wallis and Dunn test, p < 0.05) between IF1 and MF2 during the NI-cold period, and %EPT between IF1 and MF1 during the I-warm period ( Figure 5).
The PCA of the hydrological variables computed over the 90 days before benthos sampling and carried out separately for the two sampling periods, explained 65.7% of the overall variability in the datasets (see Table S4 and Figure S1). Since many variables were significantly correlated, only two parameters were selected for the RDA. Specifically, FRE LF (i.e., the number of low-flow days) considerably contributed to PC1 and was significantly correlated to most of the other variables, and Q min (i.e., the minimum flow value) considerably contributed to PC2. Both parameters characterized low flows at the study sites.
The RDA demonstrated that the variation of benthos metrics can be explained by physico-chemical parameters and by the considered hydrological variables in both sampling periods (I-warm = 63%, and NI-cold = 77%, Figure 6). A positive correlation between Q min and H was detected as well as between nutrient load (mainly NO 3 − −N) and LC and between the organic load expressed by BOD 5 and TD c . Some differences between the two sampling periods could also be observed. For instance, in the I-warm period, %EPT was negatively correlated with NO 3 − −N while in the NI-cold period it was positively correlated with BOD 5 and Q min . Overall, the outputs of the regression models were consistent with the RDA results (Table 3). Particularly, the link between Qmin and H was confirmed, and the association between a benthos metric and a hydrological variable emerged. Specifically, in the I-warm period, N was positively correlated to DURLFmax, and TDc to INCL. Moreover, in the NIcold period, S was negatively correlated to QP75, EPT and %EPT to ΔQ, and C/D to INCL, while LC and TDc were positively correlated to ΔQ and INCMax, respectively. Nevertheless, correlations between NO3 − −N and %EPT (negative) and LC (positive), between A/H and T (negative) in the I-warm period, and between pH and N (positive) in the NI-cold period were also detected. Overall, the outputs of the regression models were consistent with the RDA results (Table 3). Particularly, the link between Q min and H was confirmed, and the association between a benthos metric and a hydrological variable emerged. Specifically, in the I-warm period, N was positively correlated to DUR LFmax , and TD c to INC L . Moreover, in the NIcold period, S was negatively correlated to Q P75 , EPT and %EPT to ∆Q, and C/D to INC L , while LC and TD c were positively correlated to ∆Q and INC Max , respectively. Nevertheless, correlations between NO 3 − −N and %EPT (negative) and LC (positive), between A/H and T (negative) in the I-warm period, and between pH and N (positive) in the NI-cold period were also detected. Table 3. Significant results of the OLS regression (best model selected according to the Akaike information criterion-AIC, constraining the result to one independent variable) carried out for each benthos metric, accounting for all the considered environmental (hydrological and physicochemical) variables. The coefficient of determination adjusted (Adj. R 2 ), the F-statistic, and model significance are also reported. I-warm = Irrigation and warm period, NI-cold = Non-Irrigation and cold period.

Effects of Climate Change on Low-Flow Periods
In the last eight decades, the Lake Maggiore outflows displayed significant variation regarding seven hydrological metrics considered of ecological interest [31]. Five out of these seven parameters are associated to low-flow periods and confirm the tendency towards more severe low flows highlighted by other authors. In particular, increasing duration of low-flow periods was also detected by Meaurio et al. [43] for the Nerbioi River (Spain), and, in general, for the southern Europe watercourses by Giuntoli et al. [4]. Moreover, the magnitude of late summer flows (i.e., September flows) in the Ticino River showed a temporal decrease, regardless of the general hydrologic conditions of the year (i.e., normalized by the annual mean), in accordance to what was already found by the European Environmental Agency [44] for European watercourses since the 1960s and to what was projected for future decades for some English [3], French [5][6][7], and Belgian [8,45] rivers.
In some instances, the time patterns of specific hydrologic parameters are clearly influenced also by the operation of the Miorina Dam. For example, the frequency of low-flow pulses, at a maximum during the 1960s, decreased in the following decades. This can be attributed to the streamflow diversion for hydropower, reduced at that time during the weekends, in periods of low water storage in Lake Maggiore. A further modification of previous management of Lake Maggiore outflows is related to the introduction of the mandatory release of MF in the Ticino River, started in 2009, and revised (i.e., MF was slightly increased) in 2015. Though this new allocation would have forced a general increase of regulated outflows during summer, the decrease of average normalized flows in September seems consistent with the mentioned climate modification. However, disentangling the effects of climate change from those of changing water management on long-term flow modifications is a challenging task, frequently requiring the reconstruction of natural (i.e., unaffected by regulation) streamflow time series [46,47]. Indeed, interannual changes in flow management by dams can occur in response to operational objectives, reservoir regulation rules, and hydro-meteorological conditions [48]. Additional change of river discharge may be determined by any increase in storage capacity in the regulated catchment as well as by land use modifications, especially in mountainous areas, due to the abandonment of traditional grazing practices [49,50].

Effects of Low Flows on Benthic Macroinvertebrates
Consistently with the long-term dataset, the analysis of the 6-year streamflow time series in the river reaches downstream of water withdrawals revealed major hydrologic alteration during the warm months. River discharge in the three study sites clearly differed for the magnitude of annual minima and for mean monthly flow in September, according to the expected gradient from the site with the largest flow reduction (MF2) to that with the smallest one (IF1). Apart from the released MF values, the two MF sites differed also for the variability of discharge during low-flow periods, lower at MF1 than at MF2, and the consequent higher values of parameters related to the rate of change, due to a flatter pattern of low flows between flow peaks. These hydrological differences did not determine relevant biological differences in overall structural and functional composition between sites. However, as expected, assemblages changed in relation to the sampling period (i.e., irrigation warm vs. non-irrigation cold period). In particular, sampling period acted as a key determinant of zoobenthic taxonomical composition but not of functional composition, due to the partial substitution by different taxa occupying the same feeding niche over the different sampling periods. Chironomidae and Ephemerella (mainly gatherers and grazers) had higher densities in the NI-cold than in the I-warm period, while Baetis and Leuctra (belonging to the same feeding groups) had opposite peaks in abundance. These patterns were related mainly to the phenology of the mentioned insect taxa. In fact, the period of maximum abundance of the aquatic life stages of insects is quite predictable, particularly for monovoltine taxa, since the occurrence of the main stages, such as emergence and eggs hatching, is determined primarily by temperature and photoperiod [51]. For the most common monovoltine taxa of the Ticino River (i.e., Hydropsychidae, Leuctra, and Ephemerella), this period is in accordance with that identified by Dolédec [52] in a study of the seasonality of benthic assemblages in the Ardeche River (France). Conversely, multivoltines have the opportunity to adapt to environmental modifications by frequent generations [53,54], so that yearly variations of abundance in non-pristine rivers are probably more linked to pressures like hydrological alteration, rather than to seasonality. The highlighted differences in taxonomical composition between sampling periods include the poorer variability among samples in the I-warm, compared to the NI-cold period. Differences in taxonomical composition, as well as in variability among samples, were already described by Quadroni et al. [36], both at unperturbed sites and at sites characterized by mild hydrological impairment partly preserving the seasonal flow variability, and they were ascribed to natural seasonal dynamics in assemblages.
Sampling period also influenced the values of biological metrics, always consistently at the two sites located immediately below water withdrawals. Conversely, some metrics showed a different pattern between MF sites and the IF site. In particular, %EPT was significantly higher in the I-warm than in the NI-cold period at MF1 (and, to a lesser extent, at MF2), but not at IF1, and this is mainly attributable to the already mentioned seasonality in the abundance of the Trichoptera Hydropsychidae (i.e., the major component of EPT in the Ticino River), particularly visible at MF1. A decrease of %EPT was already reported in response to water diversion as an effect of the increase of less-sensitive non-EPT taxa [55]. However, Wills et al. [56], in a comparative study of different degrees of flow reduction, found that the density of EPT taxa in condition of mild diversion (50% flow reduction) was higher than at a reference unperturbed site, while it was lower when a flow reduction of 90% occurred. Indeed, in the Ticino River, we observed an average increase of the %EPT at MF sites compared to IF1 during the I-warm period and a decrease during the NI-cold period, when the low-flow duration was longer, while there was almost no variation at IF1. Even if to a lesser extent, the relative abundance of multivoltine organisms (i.e., LC) was lower during the NI-cold period at the two more hydrologically impacted sites, while only minor difference between sampling periods was observed at the site characterized by residual basin contribution.
Overall, diversity (i.e., H) and richness (i.e., S and EPT) metrics were higher on average during the I-warm compared to the NI-cold period and followed an increasing trend from the site with lower discharge minima to that with flow increased by the residual basin contribution. The increase in the abundance of multivoltine taxa and the decrease in Shannon-Wiener index, taxa richness, and EPT richness were already indicated among the effects of reduced discharge on benthic macroinvertebrates [36,[55][56][57][58], although the response of taxa richness to hydrologic impairment reported in the literature is not univocal as an effect of other environmental factors affecting this community metric [59,60].
Although the degree of flow reduction was often indicated as a determinant of the studied benthos metrics, the benthos impairment in the Ticino River appeared particularly evident during NI-cold, though the MFs are slightly larger than in summer. Nevertheless, in the NI-cold period, the duration of low flows is significantly higher than during the I-warm, due to naturally lower water availability.
Redundancy analysis and regressions clarified the role of the different environmental (water physico-chemical and hydrological components) facets of low-flow periods on zoobenthic assemblages. In addition to decreasing available habitat, the reduced water availability can influence water temperature and chemistry [61,62]. Accordingly, previous studies [18,63,64] demonstrated that variations in macroinvertebrate communities cannot be explained uniquely by hydrology; other predictors require comparable attention.
Although in the investigated section the water quality of the Ticino River can be considered good, i.e., well-oxygenated water without relevant nutrient and organic pollution, our results confirmed the influence of water physico-chemical parameters on the benthos metrics during the low-flow periods. Particularly, in the I-warm period, relatively high nitrate nitrogen concentrations decrease the relative abundance of EPT, favoring at the same time multivoltine taxa such as Chironomidae, i.e., a family including species characterized by relatively high tolerance to pollution [65,66]. This may be partly linked to the increasing diffuse pollution related to irrigated agriculture, typical of the warm season [18,67]. Moreover, the degree of autotrophy of the investigated system is strictly related to the summer increase in water temperature that can be considered a proxy of the seasonal succession of instream primary producers. However, our results showed also significant correlation between benthos metrics and the hydrological conditions preceding benthos sampling. First and most importantly, the key role of the released MF in shaping the community was evidenced: A higher diversity of the community expressed by the Shannon-Wiener index was explained by higher Q min (i.e., the MF). Moreover, the main hydrological determinants of zoobenthic assemblages were related to the variability/stability of low-flow periods, in terms of number of days of low flows, difference of flow magnitude, and rate of change of discharge between consecutive days. The total density was positively correlated to the duration of low flows. Similarly, in another Italian regulated lowland river, benthos sampled after periods of increasing duration of stable flows (equal to MFs) showed a strong increase in overall density [18]. Low flows can also sustain the development of richer assemblages, while flows of high magnitude can negatively affect family richness, partially in contrast to expectations [36,55]. Difference of flow magnitude controlled the relative abundance of EPT and multivoltine taxa during the NI-cold period, comparably to what was previously reported for increases of nitrate nitrogen in the I-warm period but in the opposite way. Finally, significant increase in top-down control took place when a relatively large increase in flow magnitude occurred shortly before sampling. Accordingly, Theodoropoulos et al. [68] found that predators are favored after flow peaks, due to high concentration of preys in the space between coarse grains, justified, in turn, by the high retention rates of organic matter at these microhabitats [69,70]. Although Chironomidae is a family that includes species of predators, they are generally disadvantaged by relatively high rates of flow change close to the sampling date.

Concluding Remarks
Low flows acted as key determinants in shaping structure and functions of macroinvertebrate communities in the investigated reaches of the Ticino River.
In particular, in river reaches downstream of water diversions subjected to MF release, we observed that a decrease in the lower streamflow explained assemblages with lower diversity and richness, mostly composed by tolerant taxa with short life cycles. The studied sites covered only a limited range of MF values (from 2% to 13% of the MANF), and this probably justifies the limited difference in assemblage composition between sites.
As suggested by other authors [71,72], high flows exert a significant control on benthic assemblages of regulated rivers and, also in our study, the presence of days with increased flow within long low-flow periods appeared to counterbalance some of the abovementioned effects of low flows.
Considering the tendency towards more severe low flows observed in the Ticino River in recent decades, and the projected future exacerbation reported by regional and global scale studies, a more sustainable water resource management that matches water use with environmental protection is urgently required.
In the specific case of regulated rivers outflowing from natural lakes, such as the Ticino River, streamflow management can also affect the lake ecosystem, mainly the lake shore zones [73], as changes in outflows determine fluctuation in lake water level. In this respect, a cooperation project is actually ongoing, aimed to define transboundary (Italy-Switzerland) shared rules for water and environmental management taking into account the complexity of the Lake Maggiore-Ticino River system (Interreg Project Parchi Verbano Ticino, https://progetti.interreg-italiasvizzera.eu, accessed on 15 January 2021), with possible interesting outcomes also for analogous lake-river systems.
Referring to the design of management measures capable to mitigate the impacts of low-flow periods on downstream benthic fauna, future development of the present study could include the possibility to release experimental flow peaks aimed at reducing the negative effects of extremely prolonged low flows (i.e., equal to MF) in the river, limiting at the same time drawbacks to the lake ecosystem. This does not substitute but should add to the already acknowledged need to reduce the extent of human water uses for irrigated agriculture and hydropower.
Moreover, possible modifications of insect life cycles and distributions due to climate change [74][75][76][77][78][79], and the previously reported effects of flow management and global changes on the Ticino River thermal regime [62], indicate the need to get deeper insight on the metabolic implications of long periods of stable low flows.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding authors upon reasonable request.