Habitat Suitability Curves for Freshwater Macroinvertebrates of Tropical Andean Rivers

Sustainable river management requires a thorough understanding of the response of aquatic biota to riverine microhabitat variability. The purpose of this study was to assess macroinvertebrate hydraulic-habitat suitability in Ecuadorian Andean rivers to support habitat modelling for sustainable ecosystem management. 597 macroinvertebrate samples were collected from ten sampling stations the Yanuncay River, Ecuador. Physical, chemical, hydraulic and habitat variables were measured/calculated. Froude number, Reynolds number, substrate index and algae coverage were major drivers of macroinvertebrate response, and were used to develop suitability curves for Baetodes, Andesiops, Camelobaetidius, Ecuaphlebia, Anacroneuria, Atopsyche, Simulium and Palpomyia using General Additive Models. Standardised density contours of taxa as functions of hydraulic and habitat variables were also developed. Taxonomic response was related to body structures/shapes and feeding habits. Baetodoes, Simulium, Anacroneuria and Atopsyche preferred fast flowing waters, and thus, they could be significantly affected in case of flow reduction. Similar habitat suitability curves were developed from the main river and the tributaries, possibly due to the short distance between the sampling stations. This study fills a major knowledge gap by developing macroinvertebrate habitat suitability curves for future physical habitat simulations and environmental flow assessments in the Andean region.


Introduction
In river ecosystems, environmental and biological factors act at multiple spatial and temporal scales, influencing the distribution patterns of aquatic biodiversity. In this regard, most studies on aquatic macroinvertebrates focus on large spatial and temporal scales [1][2][3], and only few on smaller scales [4,5]. Ecological or hydrological processes acting at larger scales shape smaller-scale habitats, resulting in a nested hierarchy of interconnected habitats where species become more and more specialised as habitat conditions become site-specific [6]. Moreover, aquatic communities in lotic ecosystems have irregular spatial distribution patterns [7] due to specific physical, hydrological and/or hydraulic conditions [4,8] that shape aquatic biota through adaptation processes resulting in highly specialised species [9]. Hence, the main hydraulic (flow velocity and water depth) and habitat (substrate type) characteristics and their combined influence have been assessed as major drivers of the spatial distribution of benthic communities at the microhabitat scale [10].
Most studies on the response of macroinvertebrates to hydraulic-habitat variability have been carried out in temperate zones, indicating that primary hydraulic and habitat drivers of biotic distribution are flow velocity [11][12][13], water depth [11][12][13], Froude number [13,14], Reynolds number [15], substrate type [14] and shear stress [14]. These eco-hydraulic studies have generated basic knowledge for the development of habitat models, incorporating various discharge scenarios that reflect different levels of anthropogenic impact [16,17]. However, in tropical zones, knowledge on the hydraulic response of aquatic macroinvertebrates remains limited [18][19][20].
While in other latitudes aquatic species are used to holistically assess the effect of natural or anthropic impacts on the ecological status of rivers [21,22], in Ecuador and other Andean countries, legislation on environmental flows is established mainly on empirical hydrological methods [23,24], without considering the influence of physical or hydraulic variability on aquatic biota. Moreover, in this region the alteration of aquatic and terrestrial habitats is increasing due to the construction of dams or water reservoirs of various sizes [25,26], which is likely to cause a disconnection of habitats, impacting in this way the aquatic communities [25].
Several methods have been applied to study the responses of aquatic invertebrates to physical and hydraulic-habitat variables. The biological response to hydraulic variability has been investigated by using techniques such as adjustment of polynomials [27], exponential polynomial functions [11] or local polynomial regression [28]. Other multivariate methods received greater attention in recent years such as the Generalised Additive Models (GAMs) [16,29,30]. Some other studies used large amounts of data with Boosted Regression Trees, Random Forests [28] and fuzzy-logic-based models [28,30]. GAMs are considered an appropriate method for estimating the habitat suitability of aquatic species [31] including macroinvertebrates [16,29,30] because it enables overcoming inherent pitfalls of other more traditional methods that do not assess uncertainty and do not consider interactions among hydraulic-habitat variables [16]. Also, GAMs are appropriate for developing non-linear and non-monotonic habitat suitability curves between the response and predictor variables [31]. These curves can be further included in eco-hydraulic models such as the System for Environmental Flow Analysis (SEFA) [32] or the River Hydraulic Habitat Simulation (RHYHABSIM) [33] for physical habitat modelling purposes, aiming at assessing the instream flow requirements of macroinvertebrates to maintain life-supporting freshwater capacity and native biodiversity.
Suitability curves are expensive to elaborate, and this restricts their use for carrying out individual instream flow assessments [16]. As an acceptable cost-effective practice, suitability curves developed in large rivers are normally transferred to other river types. However, assuming similar habitat suitability of aquatic organisms in different rivers might not be correct since their response to habitat variability can vary importantly in space and time based on river size and flow variation [34][35][36]. In Ecuador, hydraulic-habitat suitability curves for macroinvertebrates do not exist; as such, this study is a pioneering effort to develop habitat suitability curves for macroinvertebrates in larger and smaller rivers.
Mitigation of flow regime alteration caused by human activities and application of stream management measures require the use of predictive models to quantify the environmental flow regimes to maintain biotic and abiotic processes within river systems [27]. However, prior to the application of such predictive models, scientific knowledge of the response of aquatic organisms to hydraulic and habitat variability is necessary [37].
The purpose of this study was to assess the response of macroinvertebrates to microhabitats characterised by certain hydraulic and habitat conditions and to develop suitability curves for future physical habitat modelling in the high Andean rivers of southern Ecuador. The specific objectives were to: (a) identify the major hydraulic and habitat variables to develop habitat suitability curves; (b) develop habitat suitability curves for the most representative macroinvertebrate taxa in the Yanuncay main river and its tributaries; (c) compare the habitat suitability curves from the main river and its tributaries; and (d) produce contour curves of (standardised) density of the most representative macroinvertebrate taxa as a function of the most influential hydraulic and habitat variables. To the best of our knowledge, this is the first attempt to develop habitat suitability curves for dominant macroinvertebrate taxa in tropical Andean lotic ecosystems.

Study Area
Andean streams are characterised by highly variable discharge with sudden peak flows, highly variable topography and only two seasons, wet and drier, whose temporal bounds may vary drastically from year to year, and with respect to the average pattern [38,39]. Henceforth, particularly this seasonality is a clear distinctive difference between tropical (Andean) and temperate mountain streams, upon which there are also distinctive characteristics in terms of biotic communities. The study area is the sub-basin of the Yanuncay river (Figure 1a), which belongs to the Paute basin ( Figure 1b). The Paute basin is one of the most important basins for the domestic and industrial development of Ecuador (Figure 1c) given that approximately 40% of the nation-wide consumed electricity is hydraulically generated in it [40]. The studied sub-basin has an area of about 418.9 km 2 ; its elevation varies between 2520 and 4080 m above the sea level (a.s.l.). The slope is significantly irregular with valleys formed by steep transversal slopes. Above 3000 m a.s.l., the sub-basin presents a dominant (i.e., 90%) grass land termed as "Pajonal" (Calamagrostis sp.) [41], a typical vegetation cover of these Andean highlands (termed as "Páramo"), as well as some Quinua (Polylepis sp.) patches (i.e., 0.78%; [42]). The central sub-basin presents evergreen high montane forest patches with land dedicated to agricultural activities and cattle raising. In the lower parts significant human settlements occur. As in most of the Paute river basin, given the particular presence of the Andean mountain range, the climate of the Yanuncay sub-basin is affected by the winds of both, the Atlantic Ocean in the east and the Pacific Ocean in the west. Upon the analysis of rainfall data in the period 1997-2018, two main rainy seasons can be distinguished (i.e., January-May and October-December) with a single three-month drier season, i.e., July-September, although with the presence of rainfall, which matches prior research conclusions [38]. variables. To the best of our knowledge, this is the first attempt to develop habitat suitability curves for dominant macroinvertebrate taxa in tropical Andean lotic ecosystems.

Study Area
Andean streams are characterised by highly variable discharge with sudden peak flows, highly variable topography and only two seasons, wet and drier, whose temporal bounds may vary drastically from year to year, and with respect to the average pattern [38,39]. Henceforth, particularly this seasonality is a clear distinctive difference between tropical (Andean) and temperate mountain streams, upon which there are also distinctive characteristics in terms of biotic communities. The study area is the sub-basin of the Yanuncay river (Figure 1a), which belongs to the Paute basin ( Figure  1b). The Paute basin is one of the most important basins for the domestic and industrial development of Ecuador (Figure 1c) given that approximately 40% of the nation-wide consumed electricity is hydraulically generated in it [40]. The studied sub-basin has an area of about 418.9 km 2 ; its elevation varies between 2520 and 4080 m above the sea level (a.s.l.). The slope is significantly irregular with valleys formed by steep transversal slopes. Above 3000 m a.s.l., the sub-basin presents a dominant (i.e., 90%) grass land termed as "Pajonal" (Calamagrostis sp.) [41], a typical vegetation cover of these Andean highlands (termed as "Páramo"), as well as some Quinua (Polylepis sp.) patches (i.e., 0.78%; [42]). The central sub-basin presents evergreen high montane forest patches with land dedicated to agricultural activities and cattle raising. In the lower parts significant human settlements occur. As in most of the Paute river basin, given the particular presence of the Andean mountain range, the climate of the Yanuncay sub-basin is affected by the winds of both, the Atlantic Ocean in the east and the Pacific Ocean in the west. Upon the analysis of rainfall data in the period 1997-2018, two main rainy seasons can be distinguished (i.e., January-May and October-December) with a single threemonth drier season, i.e., July-September, although with the presence of rainfall, which matches prior research conclusions [38].

Abiotic Monitoring
Aiming at inspecting the eco-hydraulic dynamics of the studied taxa as a function of channel type, tributary and main river channels were considered separately, as well as a "global" (i.e., simultaneous) consideration of both channel types.
Ten sampling stations (five in the main river channel and five in the tributaries), surrounded by the same riverine forest, were established in the Yanuncay sub-basin (Figure 1a). The main river stations are distributed between 2919 (Y1) and 3194 (Y5) m a.s.l. Their substratum is dominated by boulder (>256 mm) and cobble (256-64 mm), and their longitudinal slope varies between 1.0% and 2.8%. The tributary stations are distributed between 2828 (T5) and 3236 (T1) m of elevation. Their substratum is dominated by boulders, cobbles and pebbles (64-4 mm) and their longitudinal slope varies between 0.25% and 6.50%.
For the 10 sampling stations, both biotic and abiotic measurements were performed in the two sampling campaigns (i.e., C1 on November 2017, average rainy period; and C2 on July 2018, average drier period). No sampling took place in periods of high peak flows owing to safety constraints and because macroinvertebrate communities are washed away during flooding. A 200 m reach segment was chosen for the monitoring of each sampling station. The choice of the segment was based on accessibility for monitoring purposes and on the representativeness of all of the mesohabitats of the study site.

Biotic Monitoring
Thirty macroinvertebrate samples were collected from different mesohabitats using a Surber net (area: 900 cm 2 ; mesh size: 250 microns). In sampling campaign C1, 30 samples were collected at 9 sampling stations while at sampling station T5 it was only possible to collect 27 samples during the rainy period when too high discharges occurred (i.e., 297 samples were collected). In campaign C2, 30 samples were collected at every one of the 10 sampling stations (i.e., 300 samples were collected). Hence, a total of 597 macroinvertebrate samples were taken in both campaigns, 300 at the main river and 297 at the tributaries. The coverage of algae (Algae, %) and bryophytes (Bryo, %) per sampling quadrant (900 cm 2 ) were visually estimated (by the same person), previously to each macroinvertebrate sampling. Macroinvertebrate samples were fixed in the field with 4% formaldehyde solution. In laboratory, all of the organism were sorted out and identified to genus level using taxonomic keys [45]. Higher taxonomic divisions were used for organisms belonging to the non-insect class, the subfamilies of Chironomidae and undistinguishable larvae of the Xiphocentronidae family.
Once all of the macroinvertebrate organisms were removed from the sample, the remaining sediment was used to determine the organic matter content in the sediment, according to the procedure of Steinman et al. [46].

Data Processing
Rare macroinvertebrate taxonomic groups, characterised by a relative abundance lower than 0.01% [47], were discarded from the statistical analysis. Furthermore, using the absolute abundance of every taxon, the density metric (ind. m −2 ) was calculated. Density was further standardised because biological data were collected at different discharges and times. Standardisation was achieved by dividing density estimates by their arithmetic mean calculated for every sampling campaign (C1 and C2) and for each channel type (main river and tributary). Hence, four groups of "standardised by mean" density (δ) data were generated, resulting from the consideration of two sampling campaigns and two channel types.
Hydraulic measurements were processed to calculate some hydraulic variables of interest. Hence, the Froude number (Fr), the ratio between the inertial and gravitational forces, was calculated as a function of water velocity, gravity acceleration and the hydraulic depth [48]. Also the Reynolds number (Re), the ratio between inertial and viscous forces, was calculated as a function of water velocity, the hydraulic radius and the kinematic viscosity [49]. Re values lower than 500 were assumed to characterise laminar flows, whilst Re > 12,000 characterised turbulent flows and 500 ≤ Re ≤ 12,000 characterised transitional flows.
Moreover, the roughness Reynolds number, an empirical variable which relates the roughness and flow parameters at which transition occurs at the surface of the roughness element of a river bed [50] was calculated as a function of the mean water velocity in the column of water, the height of the surface roughness element and the kinematic viscosity [50,51]. Roughness Reynolds number was used in this study after [52] detecting that this empirical criterion was importantly related to the spatial variation of invertebrates. Since roughness Reynolds number describes the roughness of the near-bed flow environment, values lower than 5 were assumed to characterise hydraulically smooth flows, whilst roughness Reynolds number > 70 characterised hydraulically rough flows (i.e., in rivers with coarse bed material) and 5 ≤ roughness Reynolds number ≤ 70 characterised transitional flows [51,52].
The percent coverage of each substrate category were used to calculate the Shannon-Wiener substrate diversity (SuD, [53]). SuD indicates the heterogeneity of the substrate. Further, the substrate index (SI) was calculated summing weighted substrate percentages [54], as follows: (0.08) % bedrock + (0.07) % boulder + (0.06) % cobble + (0.05) % pebble + (0.04) % granule gravel + (0.03) % sand. Given the summation properties of its mathematical expression SI varies between 3 (for a 100% sand content) and 8 (for a 100% bedrock content) and may adopt any real value in that range of variation. Further, the larger SI, the larger the sampled substrate sizes are.

Statistical Analysis
Spearman's correlation analysis was performed on the 21 hydraulic and habitat variables for the datasets of the main and tributary river reaches. Three variables in the case of the main reach (i.e., velocity measured at 60% of the water depth, the roughness Reynolds number and the boulder substrate material) and 2 variables in the case of tributaries (i.e., velocity measured at 60% of the water depth and the roughness Reynolds number), presenting a significant correlation (r > 0.80, p < 0.05) with other hydraulic/habitat variables of higher interest for this study, were discarded from the rest of the analysis. Further, Bryophytes coverage was also discarded since it was not present in the main river reach. The Detrended Correspondence Analysis (DCA) was applied for selecting the appropriate multivariate model that was used, in turn, for determining the hydraulic and habitat variables influencing macroinvertebrate communities at microhabitat scale. DCA, based on the length of gradient (<3), suggested the use of the linear response multivariate model (Redundancy Analysis, RDA; [55]). Two RDA multivariate models were constructed [56]; one for the main river and another one for tributaries. The multivariate (DCA and RDA) analyses were carried out using non-standardised taxa densities, which, previously, were transformed using the log (x + 1) expression and standardised by species [55]. Forward selection was employed to select the hydraulic and habitat variables that had influence on the macroinvertebrate communities. For a better visualization of the resulting plot, the number of taxa was restricted to the 15 that had a best fit in the RDA models. These statistical analyses were carried out with the help of the CANOCO software version 5.02 [57].

Development of Habitat Suitability Curves
All variables selected by the RDA were used as predictors in GAMs but only the results of four hydraulic and habitat variables (Fr, Re, SI, and Algae) are presented herein. They were selected on the basis of their significance in the RDA, their similar use in past research [15,16], and the more meaningful interpretation of the respective GAMs based habitat suitability curves. Eight taxa were selected as response variables in the GAMs analysis based on their high occurrence and dominance in the samples, their sensitivity to river water quality [58], their significance in the RDA and their meaningful response in the GAMs. The application of the GAMs was achieved through the use of the Habitat Selection program (version 3.0, Jowett Consulting Ltd., Pukekohe, New Zealand) with a natural logarithmic (ln) link function [16,59]. In this way, the general expression of a full model is as follows: where Y is the dependable variable (standardised density, δ, in the current case); X i is the i-th hydraulic/habitat variable, with i being an integer sub- A test for non-linearity was carried out to check on whether non-linear GAMs are needed to depict the biotic responses to hydraulic/habitat variables. An effective degrees of freedom value, df = 3, was used. This value seemed to produce a realistic degree of smoothing of the resulting curves and followed reasonably well the average evolution of the plotted observations. A Poisson-type error distribution was considered, potentiated by a further over-dispersion correction [16], which occurs when the sample's variance is greater than that of a Poisson distribution (i.e., σ 2 = 1). The over-dispersion correction prevents standard errors from being underestimated and significance of model variables being overstated [16,59].
To measure the accuracy of the resulting GAMs, the goodness of fit index "total deviance explained" [16,59] was used. In this regard, the residual deviance is a measure of discrepancy between the observations and the GAMs fitted values. It is computed using the likelihood of the GAM of interest as a function of the adopted error distribution. The null deviance on the other hand is the deviance for a model with just a constant term. The total deviance explained, expressed either as a proportion or as a percentage, is the relative difference between the residual and the null deviances. Henceforth, it is calculated as the quotient among the difference between the null and the residual deviances, divided by the null deviance.
For depicting the results of this analysis, plots were prepared combining the observations of the dependable variable (δ), linked to the primary Y-axis, with the "Contribution" (expressed in ln(δ) units) of the X i hydraulic/habitat variable to the full GAM, linked to the secondary Y-axis. Further, standardised density contours of the different taxa as a function of hydraulic and habitat variables, concretely, Fr, Re, SI and Algae, were also developed. In every contour plot, the ranges of variation of two of these variables were inspected while fixing the values of the other two remaining variables to their "optimal" value.

Environmental Characteristics and Macroinvertebrate Taxa
Water temperature (between 8.2 and 14.0 • C) and pH (between 6.6 and 8.1) were relatively similar at all sampling stations (Table 1). Average turbidity (<8.5 NTU) and total dissolved solids (<0.1 g·L −1 ) values were low. Dissolved oxygen values were relatively similar in the main river and tributaries (6.6-10.5 mg·L −1 ). Conductivity was slightly lower in the main river (<95.3 µS·cm −1 ) than in the tributaries (<129.0 µS·cm −1 ). Water discharge values varied between 0.78-3.32 in the main river and 0.02-1.39 m 3 ·s −1 in the tributaries. Similar ranges of variation of the following variables were observed/estimated at/for the main river and the tributaries: 3.0 ≤ water depth ≤ 45 cm, 0.0 ≤ water velocity ≤ 205 cm·s −1 , 0.0 ≤ Fr ≤ 2.3 and 0.0 ≤ SuD ≤ 2.0.
In both channel types, the dominant substrate types were boulders (average in main river = 26.9%; average in tributaries = 21.0%), fine cobbles (average in main river = 15.2%; average in tributaries = 17.0%) and very coarse pebbles (average in main river = 14.9%; average in tributaries = 17.4%). Re was higher in the main river (average Re = 107,204) than in the tributaries (average Re = 75,833). This was also the case of roughness Reynolds number, as the main river (average) roughness Reynolds number value was approximately 2.38 × 10 6 , whilst the tributaries (average) roughness Reynolds number value was about 2.08 × 10 6 , which is in agreement with the average characteristics of the substrate material of both channel types. Average organic matter content in the sediment was lower in the main river (9.9 g·m −2 ) than in the tributaries (37.9 g·m −2 ). Average Algae were similar in the two channel types (main river = 20.3%; tributaries = 16.1%).
In total, 61 macroinvertebrate taxa were recorded in the main river and tributaries ( Table 2). Limonicola was the only taxa that was present only in one river type (i.e., main river). The most abundant taxa were Chironomidae, Oligochaeta, Baetodes, Metrichia and Andesiops contributing with 91.7% to the total relative abundance and 30 taxa had less than 0.002% of relative abundance.

Variables Influencing the Macroinvertebrate Communities at Microhabitat Scale
The RDA identified five variables (Fr, SuD, Algae, organic matter content in the sediment, Re) that affected the spatial distribution of the macroinvertebrate communities in the main river ( Figure 2; Table 3). The first axis of RDA explained 10.4% of the variation in the data set and showed a gradient related to Fr and Re. The second axis explained 3.9% of the variation and indicated a gradient related to Algae (Figure 2a; Table 3a). RDA selected 7 significant hydraulic and habitat variables (Fr, fine cobble, SI, Re, SuD, Algae, Bryo) in the tributary dataset. The first axis of RDA explained 6.6% of the variation in the data set and reflected a gradient related to Fr, Re, SI and Bryo.
Water 2020, 12, x FOR PEER REVIEW 9 of 21     The second axis explained 4.1% of the variation and indicated a gradient related to fine cobble and SuD (Figure 2b; Table 3b). Palpomyia, Molophilus, Austrolimnius, Oligochaeta taxa were negatively associated with faster flow habitats (i.e., characterised by higher Fr values) in both, the main river and tributaries. Further, in the tributaries they are negatively associated with SI. In contrast, Simulium and Baetodes taxa were positively associated with Re in both channel types. Ecuaphlebia and Hydracarina taxa were positively related to SuD in the tributaries while Hydracarina and Leptohyphes were positively associated to SuD and organic matter content in the sediment in the main river. Leptohyphes, Andesiops, Atopsyche and Anacroneuria were positively related to fine cobble substrate in the tributaries. Limnophora, was positively correlated to Algae in tributaries; similarly, Phylloicus is positively correlated to Algae in the main river.

Individual Taxa Response to Selected Hydraulic and Habitat Variables
The application of the GAMs method on the standardised densities (δ) of the selected 8 taxa in the Yanuncay River (i.e., global analysis) using only significant variables explained different percentages of the total deviance (Table 4). In the case of Baetodes and Simulium, the explained percentages of the total deviance were the highest, 34.5% and 32.4%, respectively, for the full GAMs; whilst in case of Ecuaphlebia and Atopsyche the four significant variables explained the lowest percentages of the total deviance (i.e., 14.7% and 13.0%, respectively). Fr explained the highest deviance in the δ of Baetodes, Simulium, Andesiops, Camelobaetidius, Atopsyche and Palpomyia in the GAM model (18.6%, 16.7%, 11.7%, 9.2%, 5.9% and 13.5%, respectively). SI explained the highest deviance in the δ of Anacroneuria and Ecuaphlebia with 8.3% and 5.8%, respectively. Further, the F test (F ratio > 2.5 and p < 0.05) suggested significant non-linear relationships between the δ of all of the 8 studied taxa and Fr, Re, SI and Algae supporting as such the use of non-linear GAMs in the current study.
Baetodes (Figure 3a (Figure 4p) did not show clear trend as response to this variable; moreover, response of these taxa to Algae varied importantly between the main river and the tributaries.
Seeking clarity, Figures 3 and 4 depict only the confidence band encompassing the global habitat suitability curves. In general, there was less confidence in the regions of the habitat suitability curves where standardised density data were scarce, indicated by increased width of the confidence bands (Figures 3 and 4). Table 2. Occurrence of taxa in the main river (Y) and tributaries (T) within the Yanuncay sub-basin, and relative abundance calculated for the global (i.e., total) data set.  Table 3. Results of the Redundancy Analysis (RDA) applied on the data of (a) main river; and (b) tributaries. Number of samples considered in the analysis, N = 300 (main river) and N = 297 (tributaries). p = significance value.  Figure 5 illustrates the standardised density (δ) contour curves of the different taxa as a function of Fr-SI and of Fr-Algae, as described by the full GAMs. High response δ of the genera Baetodes (Figure 5a), Simulium ( Figure 5b) and Anacroneuria (Figure 5c) occurred in habitats characterised by Fr > 0.5 and larger substrate sizes. Atopsyche (Figure 5g) was present in habitats with Fr > 0.5 but smaller substrate sizes. Andesiops (Figure 5d) Camelobaetidius (Figure 5e), Ecuaphlebia (Figure 5f) happened in habitats with Fr < 0.5 and larger substrate sizes. The Palpomyia genus (Figure 5h) occurred in quiet habitats (Fr < 0.25) and smaller substrate sizes.

Discussion
The current study analysed the response of aquatic macroinvertebrates to hydraulic and habitat variability occurring at the microhabitat scale. Macroinvertebrate communities were sampled at habitats distributed along the main and tributary channels of the Yanuncay River, at different flow conditions. 597 samples were included to obtain robust results. Environmental conditions at the macro-and micro-scales were measured/calculated. At the macro-scale, physical-chemical variables, slope, etc. can represent filters for the composition of macroinvertebrate community [10] and for the distribution of certain taxa; similarly to the role of habitat and hydraulic variables at the microhabitat scale.

Variables Influencing Macroinvertebrate Distribution at the Microhabitat Scale
Hydraulic and habitat characteristics (Fr, Re and substrate size) are important for determining the spatial distribution of aquatic taxa at the micro-scale [14,15,54]. The distribution of taxa in rivers depends on the capacity of its organisms to occupy different habitats, which, in turn, depends on their body structure (i.e., nails, gills, suction cups, etc.) or shape (hydrodynamic, cylindrical or flattened) developed to persist in a particular habitat [60][61][62].
In this study, Simulium, Anacroneuria and Baetodes genera were found in habitats with faster and turbulent flows, which is likely the result of their specific adaptations to this environment. The

Discussion
The current study analysed the response of aquatic macroinvertebrates to hydraulic and habitat variability occurring at the microhabitat scale. Macroinvertebrate communities were sampled at habitats distributed along the main and tributary channels of the Yanuncay River, at different flow conditions. 597 samples were included to obtain robust results. Environmental conditions at the macroand micro-scales were measured/calculated. At the macro-scale, physical-chemical variables, slope, etc. can represent filters for the composition of macroinvertebrate community [10] and for the distribution of certain taxa; similarly to the role of habitat and hydraulic variables at the microhabitat scale.

Variables Influencing Macroinvertebrate Distribution at the Microhabitat Scale
Hydraulic and habitat characteristics (Fr, Re and substrate size) are important for determining the spatial distribution of aquatic taxa at the micro-scale [14,15,54]. The distribution of taxa in rivers depends on the capacity of its organisms to occupy different habitats, which, in turn, depends on their body structure (i.e., nails, gills, suction cups, etc.) or shape (hydrodynamic, cylindrical or flattened) developed to persist in a particular habitat [60][61][62].
In this study, Simulium, Anacroneuria and Baetodes genera were found in habitats with faster and turbulent flows, which is likely the result of their specific adaptations to this environment. The resistance capacity of Simulium to faster flows is related to the presence of a circular hook structure in the last body segment [63], which allows the organism to be fixed to the bottom of the river bed. Likewise, organisms of the Baetodes genus have hydrodynamic shape, small size and tarsal nails to hold onto the substrate, allowing significant resistance to faster currents [64]. Similarly, the organisms in the genus Anacroneuria have strong and flattened body structure [65] against strong current, while Camelobaetidius has a hydrodynamic shape and cushion legs [66]. Both, Anacroneuria and Camelobaetidius, typically occur in habitats with fast currents; nevertheless, in this study Camelobaetidius preferred slower flows.
Palpomyia, Molophilus, Phylloicus and Oligochaeta are characterised by the absence of specific body structures or adaptive body form to high currents; hence, they are prone to be negatively affected by faster flows. Further, these groups are frequent in low flow areas (pools) where fine sediment and OM accumulate [7,10,67], which was also observed in this study. Herein, Huleechius were present in habitats of the main river with slower flow and higher OM, in contradiction to the results from the study of Ríos-Pulgarín et al. [64] applied on a high mountain river. Furthermore, Austrolimnius was also present in the main river and tributary habitats characterised by lower flow and higher OM conditions, which is similar to what was found by González-Trujillo and Donato-Rondon [68] that carried out a study above 3500 m a.s.l.
Another influential variable in the habitat distribution of certain aquatic taxa is the diversity or heterogeneity of substrate (SuD). Generally, organisms prefer substrates of different sizes, which provide them better refuge against different flow conditions or predators [69,70]. RDA showed that Hydracarina was associated to these refuge habitats with diverse substrate from which organisms could find further refuge in the hyporheic zones at extreme (high or low) flow events [71]. Ecuaphlebia, belonging to the Leptophlebidae family, was also present in heterogenic substrate likely using these sites as refuge, similarly to what was observed in other latitudes by Dostine et al. [71] with respect to the Leptophlebidae family. On the other hand, Cobb et al. [72] refers to the fact that large substrates provide stable habitat to aquatic insects; this has been observed in this study specially in the tributaries for Leptohyphes, Andesiops and Atopsyche that preferred the (128-64 mm) thick substrate.

Individual Taxa Response to Selected Hydraulic and Habitat Variables
The above refers to the RDA results that outlined linear responses of certain taxa to hydraulic and substrate variables. Furthermore, habitat suitability curves of taxa for certain microhabitat conditions were developed using GAMs. Herein, the observed preference of the Simulium genus for sites with faster and turbulent flow is consistent with the results of the study of Palmer et al. [73]. Simulium fixes its body to a thick substrate and takes advantage of faster flow by filtering water for capturing suspended food particles; however, in small channels these organisms are opportunistic and quickly adapt to different conditions [74]. Anacroneuria is also well adapted to faster flow and turbulent habitats as indicated by the respective habitat suitability curve. In Andean rivers Anacroneuria has feeding preference for the genus Simulium [75]. However, in this study Anacroneuria standard density decreased more rapidly than the Simulium standardised density under faster flow conditions; hence, these habitats could represent refuge for this genus against predation.
Baetodes are herbivorous [76]; they were present in habitats with faster flow (i.e., Fr > 0.5), which results in a lower coverage of filamentous algae. This genus is likely to have a feeding habit specialised on algae species resistant to fast flow. Atopsyche is mainly predator [77] and its feeding strategy is supported by faster flow but less turbulent habitats with large substrate sizes. Ecuaphlebia genus was present mainly in lower velocity habitats with intermediate substrate size, which contradicts the findings of Vimos-Lojano et al. [78] that refers to Ecuaphlebia being present in habitats with faster currents and larger substrate in an Andean subcatchment with elevations higher than 3500 m of elevation. The genus Andesiops, despite having a hydrodynamic form for swimming and tarsal nails for attaching to the substrate [77] do not frequent faster and turbulent flow environments. Population-genetic structure of Andesiops was studied by Finn et al. [79]; however, more detailed studies of this genus are needed, given the dominance of this genus in the Andean region and its potential ecological importance.
The low resistance to high flows of Camelobaetidius and Palpomyia is mainly due to their biological features that are poorly adapted to strong currents [77,80]; hence, these groups are frequent in areas of pools or slow flow. According to Tomanova et al. [77] Camelobaetidius can be scraper or collector-gatherer; in the Yanuncay sub-basin this genus most probably has the latter behaviour because it is present in habitats with low algae coverage. Palpomyia could be either predator, scraper or collector-gatherer [77]. In this study, this genus occupied habitats with small substrate and intermediate or large algae coverage, which suggests scraper and collector-gatherer feeding habits.
The current study revealed very few differences in hydraulic and habitat responses of the studied taxa for Fr, SI and Re as a function of channel types (i.e., main river and tributaries). These findings contradict other studies. Jowett [34] indicated a scaling effect in hydraulic responses for macroinvertebrates in New Zealand rivers, with preferred depths and mean velocities in small streams being considerably lower than those noticed for large rivers. Furthermore, Mérigoux et al. [81] revealed differences in hydraulic responses of macroinvertebrates between small and large streams; as such, they questioned the transferability of models of hydraulic habitat responses between rivers. In the current study, the relatively short spatial distance between the sampling stations located at the main river and tributaries and their similar characteristics are likely to have resulted in studied taxa exhibiting very similar responses to the previously mentioned hydraulic and habitat variables.
Hydropower projects are rapidly increasing in the Neotropics. In the Andean region, around 50% of new dams are classified as high impact (i.e., major break in river connectivity, deforestation, threaten to biodiversity, etc.) and just 19% as low impact [25,82]. Owing to climate change, glaciers in the tropical Andes have been retreating for the past several decades, leading to a shortage of water supply downstream [83]. Moreover, in the Paute river basin, where the study site is located, increasing peak flows and decreasing low flows are expected in the future because of climate change [84].
A dam construction is planned in the near future in the upper part of the study sub-basin as well as climate and land use changes are expected to happen in the mid-term, which will likely result in hydrological/hydraulic alteration in the studied river with various consequences on associated ecosystems. Henceforth, macroinvertebrate taxa that need higher discharge, faster flow and more turbulent flow conditions might disappear (i.e., Baetodes, Simulium, Anacroneuria and Atopsyche) to give space to taxa that are linked to slow flow conditions (i.e., Palpomyia, Oligochaeta, Molophilus).
Further, since flow regime controls the occurrence of algae, dam regulated flow could lead to an increase of biomass and density of some filamentous algae [85] and, as such, might impact macroinvertebrate communities, for instance, increasing the abundance of taxa, which have scraper feeding habit, or might result in the change of feeding habits. For instance, Camelobaetidius is most probably collector-gatherer in the studied river but with different algal coverage it might shift its feeding habit to scraper. Algal coverage change might alter ecosystem functioning by modifying the distribution and abundance of the base of the aquatic food web [86].
To maintain hydrological/hydraulic conditions for protecting river biota and good ecological status of the studied river, direct physical changes including sand extraction or river straightening need to be strongly under control. Further, in the view of future hydrological/hydraulic changes owing to dam construction and operation, environmental flow regime needs to be established through the incorporation of data on flow-ecology relationships into operational use of environmental flow definition [85].

Conclusions
This study was carried out in the Yanuncay sub-basin, southern Ecuador. However, despite the fact that results achieved cannot be directly transferred, they may be considered in other similar studies applied in Andean aquatic systems. Knowledge of the hydraulic-habitat response of macroinvertebrates at the microhabitat scale is important for a more holistic determination of environmental flows; however, this knowledge is still poor in the Andean context.
Eight hydraulic and habitat variables (SI, Fr, Re, Algae, fine cobble, SuD, Bryo and organic matter content in the sediment) were regarded as influencing the distribution of aquatic macroinvertebrate communities in the main river and tributaries. Further, GAMs were used to depict microhabitat responses of several macroinvertebrate taxa in terms of these hydraulic and habitat variables. The analysis showed more meaningful habitat suitability curves only for Fr, Re, SI and Algae. Eight taxa were selected as response variables in the GAMs analysis based on their high occurrence and dominance in the samples, their sensitivity to river water quality, their significance in the RDA and their meaningful response in the GAMs. Baetodes and Simulium had flow preferences represented by Fr > 1.0 and larger substrate size while Anacroneuria had similar habitat preference for the substrate though for slightly slower flow (0.5 < Fr < 1.0). Andesiops, Camelobaetidius and Ecuaphlebia preferred slower flow (0.5 < Fr) though still large substrate. On the other hand, Palpomyia had affinity to lentic areas with small substrate size. Algae influenced the presence of certain taxa (i.e., Baetodes, Camelobaetidius, Palpomyia) in connection with their scraper or collector-gatherer feeding habits. Atopsyche is mainly predator and its feeding strategy is supported by faster flow but less turbulent habitats with large substrate sizes. Baetodes, Simulium, Anacroneuria and Atopsyche taxa could be affected in case of discharge reduction by human activities (i.e., the presence and operation of dams). Palpomyia could be affected by sudden increase of discharge resulting from natural flood events or water released from dams.
The results of this study contributed to the scientific knowledge of the hydraulic responses of aquatic organisms, which is necessary to design predictive models for quantifying the necessary flow regimes aiming at maintaining biotic and abiotic processes within river systems. Future studies should focus on physical habitat modelling and environmental flow analysis using the eco-hydraulic knowledge of aquatic taxa developed in this and similar studies, as a contribution to better water resources management in the Andean region. In general, there was less confidence in the regions of the generated habitat suitability curves where standardised density data were scarce. This should be considered when using the generated habitat suitability curves in the scope of further analyses, particularly, habitat modelling. Funding: This research was funded by the Dirección de Investigación de la Universidad de Cuenca (Research Directorate of the Universidad de Cuenca, DIUC), grant for the project "Efectos del estrés hídrico sobre la biodiversidad en ríos abastecedores de la ciudad de Cuenca (Ecuador)". The APC was kindly waived by MDPI.