Taxonomic and non-taxonomic responses of benthic macroinvertebrates to metal toxicity in tropical reservoirs. The case of Cantareira Complex, São Paulo, Brazil

: Benthic macroinvertebrates are organisms that are recognized as water quality bio-indicators. A wide variety of indices and metrics have been shown to respond to a variety of anthropogenic impacts, usually under a general condition of environmental impairment. The absence of a clear distinction in the relations between specific pollutants and biotic variables is very common and can lead to biased interpretation of biomonitoring. The aims of this research were to test taxonomic and non-taxonomic responses to specific environmental conditions instead to general conditions. For this purpose, we estimated the theoretical toxicity by comparing toxicity values published by EPA with metal concentrations in water and sediments. Then we tested the responses of biological variables to toxicity and other environmental conditions using the linear mixed effects models approach. We generated 32 models considering 24 different biological metrics and indices that were grouped in five levels. Taxonomic and abundance metrics were best predictor than functional or tolerance-based indexes. The strongest model was that which considered subfamily taxonomic resolution responding to Al_w and Cr_s.


INTRODUCTION
The assemblage of benthic macroinvertebrates (BMI) is made up of heterotrophic macroscopic organisms, which have a direct relation with sediments. In general, these organisms have low mobility and are primarily detritivores (Carvalho & Uieda 2009). BMI are fundamental to the maintenance of aquatic ecosystems because they recover matter from sediment debris (Frauendorf et al. 2013, Whiting et al. 2011. BMI promote sediment oxygenation through bioturbation, which in turn affects biogeochemical cycles (Puigagut et al. 2014, Shang et al. 2014).
These organisms are widely used as bioindicators in continental waters (Czerniawska-Kusza 2005, Zilli et al. 2008, Azevêdo et al. 2015, Erba et al. 2015, Damanik-Ambarita et al. 2016 because they have many characteristics that are desirable in a bio-indicator. BMI are widely distributed throughout virtually every aquatic ecosystem (Bonada et al. 2006). Their low mobility means that BMI cannot escape from any damaging agents, and this characteristic permits us to identify environmental gradients in pollution (Beghelli et al. 2014) or other gradients of change, for example (Brand & Miserendino 2015). Furthermore, BMI rapidly respond to environmental alterations, and they include a wide variety of organisms with distinct ranges of tolerance to pollutants (Mandaville 2002, Bonada et al. 2006. BMI´s relations to sediments lead them to respond to problems that are not usually detectable in water because pollutants tend to accumulate on the bottom over time, especially in lentic systems (Bettinetti et al. 2012). Because of the position of BMI in the food webs, these organisms receive matter from other trophic levels such as phytoplankton (Yu et al. 2013), macrophytes (Silveira et al. 2013), zooplankton (Tang et al. 2014) and riparian vegetation (Scharnweber et al. 2014). Given their position, any stress that occurs in the water, surrounding land or sediments has the potential to cause a BMI response.
Potentially toxic metals are one of the most worrisome pollutants in continental waters (Zhang & Liu 2014, Schneider et al. 2014, López-Dovál et al. 2015. Attention to these metals is necessary, primarily because they are widely distributed, they are persistent in the environment and they can be bioaccumulated, leading to damage in ecosystems as also to human health (Schwarzenbach et al. 2006, Bin et al. 2013, Alves et al. 2014, Roig et al. 2016. Furthermore, analytical techniques for determining metals in waters (US-EPA 1992) and sediments (US-EPA 1996) are well-established, and metal toxicity levels are known (US-EPA 1999).
In considering BMI assemblages, there is a vast variety of metrics, indices and taxa that are described as bio-indicators in the literature. Among these measurements, there are classes of indices, namely, the index directly related to diversity (Azevêdo et al. 2015), the proportion of tolerant or sensitive organisms within the assemblage (Beghelli et al. 2012, Clews et al. 2014, the BMI density of organisms, and indices related to functional groups (He et al. 2015).
In addition, the majority of studies that have focused on BMI as bio-indicators were designed to distinguish between damaged and non-damaged areas to relate the general environmental quality to a variety of environmental variables. These studies primarily relate "damaged" or "non-damaged" conditions to the given land use, landscape ecology, organic pollution or general toxicity (Clews et al. 2014, Fanny et al. 2013, Gonçalves & Menezes 2011). The situation is no different in Brazil, where the use of BMI as bio-indicators is focused within lotic systems (Baptista 2008, Baptista et al. 2007, Couceiro et al. 2012, Egler et al. 2012, Ferreira et al. 2011, Mugnai et al. 2008, Junqueira & Campos 1998, Suriano et al. 2011. Despite the known relevance of reservoirs in Brazil (Tundisi & Matsumura-Tundisi 2008), there is a clear lack of information about the use of benthic macroinvertebrates as bio-indicators in these environments. As a contribution, the Companhia Ambiental do Estado de São Paulo, or CETESB, developed a multimetric index for littoral benthic macroinvertebrates for reservoirs. This agency tracks the taxa richness, the Shannon-Wiener index, the dominance of tolerant groups, the number of sensitive taxa and an index of sequential comparison (CETESB 2012).
These studies have made a great contribution to the fields of Ecology and Ecotoxicology, especially in view of the fact that ecosystems are complex and that the biomonitoring of general conditions is necessary. However, approaches that consider specific classes and attempts to isolate the specific effects of distinct pollutants or damaging agents are also desirable to increase our knowledge about their effects in nature and also to better delimit the bio-indicator potential of different organisms. Once a general environmental status approach is chosen, it is difficult to discriminate among the damaging agents to which the BMI are really responding or if there are any synergistic or competing effects (Zhang & Liu 2014).
An important contribution was made by Zhang & Liu (2014), which provided an ecological risk approach to metals through BMI metrics (functional and tolerance measurements). These authors assessed the effects of metal toxicity on benthic macroinvertebrates in Lake Baiyangdian, China along a gradient of human impacts. Studies that take this direction are desirable because they provide an important complement to the more general approach, and because they allow for predictions that are more accurate. Furthermore, distinct analyses of specific climatic regions are also complementary and beneficial.
By testing BMI responses to metal toxicity, taxonomic and non-taxonomic metrics and indices are expected to vary according to the metal toxicity in water and sediments (Zhang & Liu 2014). Based on this assumption, the present work is premised on the following hypotheses: differential habits and sensibilities may promote distinct responses. Metrics and indices related to the proportion of tolerant groups (i.e., dominance indices and indices based on tolerant taxa) may generally be favored by high toxicity. The taxa density and richness are expected to decrease with increasing toxicity, and a finer taxonomic resolution is expected to be more sensitive than the coarser measurements. Metals that present high theoretical toxicities are also predicted to have a wider influence on BMI assemblages, and various indices and/or metrics can respond distinctly to the toxicity exerted by different elements.
In considering niche theory and its relation to interspecific competition, higher toxicities must promote a decrease in the relative abundance of sensitive groups and an increase in the most tolerant ones. As a result, complementary relations may be observed, with increases in the relative abundance of one group coupled with a decrease in another as a response to the same toxicities (Colas et al. 2014).
The aim of this study was to test the mixed effects of metal toxicity in sediments and water over several BMI assemblages and their responses according to indices and metrics based on distinct levels in the taxonomic resolution, functional indices, abundance and tolerance. With this outlook, we aim to provide models that will generate parameters for use in biomonitoring or future research.

Study area
The Cantareira Complex is formed by five reservoirs, namely Jaguari, Jacareí, Cachoeira, Atibaia and Paiva Castro, which are linked by tunnels and channels. For practical purposes, Jaguari and Jacareí can be considered as the same reservoir once there is no physical barrier between them. Each reservoir receives water from the former reservoir and from a primary tributary river ( Figure 1). Waters from the Cantareira Complex are used to supply the people in the Metropolitan Region of São Paulo (RMSP), Brazil. This complex provides approximately 30% of the water consumed by 20 million people in the RMSP according to data published by the company that is responsible for water treatment and supply, which is known as the SABESP (SABESP 2015).
The Cantareira Complex is surrounded by an urbanized area that was originally occupied by Atlantic rainy forest. Based on satellite imagery, Whately & Cunha (2007) concluded that only 21% of the original forest cover was preserved until 2003. The majority of the Cantareira Complex surroundings was modified into agriculture fields, pasture and abandoned fields. With respect to the area of each reservoir's watershed, the highest proportion of natural vegetation (37%) occurs in the Paiva Castro watershed.
The Jaguari reservoir watershed has 17.7% natural forest coverage, and Cachoeira has 17.5; Jacareí covers just 11.6%. Furthermore, regional industrialization, mining activities and the presence of extensive roads are potential sources of damage that may contribute to water quality degradation (Whately & Cunha 2007 (Valderrama 1981) in the Cantareira Complex reservoirs in May and June 2013, at sampling points next to those chosen for the present work. In considering the proximity with the sampling points from the present study, Figure 2 shows the approximate concentrations of total phosphorus and chlorophyll-a determined by this study, and also the trophic state index by Carlson as adapted by Lamparelli 2004 for tropical reservoirs. These data indicate the oligotrophic state of all of the sampling points during the period. Higher trophic conditions occurred in the Jaguari-Jacareí reservoir.
In relation to the sediment, a recent report by the Companhia Ambiental do Estado de São Paulo, or CETESB for 2014 indicated that acute toxicity occurred in Vibrio fischeri in more than 80% of the sediment samples from those reservoirs. The sediment toxicity was considered high in more than 50% of the samples analyzed by CETESB 2015.

Sampling points
Samples were collected during July and August, 2013 at three points from a depth of 3 m (3.23 ± 0.5) ( Table I). The sampling points were distributed along reservoirs that were divided into approximate head, middle and dam regions.

Benthic Macroinvertebrates (BMI)
Sediments were collected for BMI counting and identification at each sampling point. The BMI sampling was performed by using an Ekman dredge (area=0.021 m 2 ) that was launched three times to make one cumulative sample. The sediment was passed through a 0.2 mm sieve and preserved in a 4% formaldehyde solution. Afterwards, the organisms were identified to the genus or species level by using identification keys and manuals (Brinkhurst & Marchese 1992, Mansur et al. 2004, Marchese 2009, Borkent & Spinelli 2007, Trivinho-Strixino 2011. After the taxonomic identifications, taxonomic and non-taxonomic metrics and indices that considered the taxonomic species/ genus, subfamily and family levels, as well as the tolerance, functional and abundance levels were calculated. The multimetric index for sublittoral BMI from CETESB 2012 was calculated as a general regional index specific to reservoirs. Only taxonomic groups that occurred in at least one sample from each reservoir were considered. The biological metrics and indices were organized by category depending on the category in use (taxonomic or non-taxonomic) and by level, depending on the given taxonomic level (Table II).

Sediment
Three samples (see Table I) from each reservoir were collected (total = nine samples) with an Ekman dredge, and organic and inorganic fractions of coarse (> 212 µm) and fine sediments were calculated by sieving the dried material. The organic matter content was determined for each fraction by calcination in a muffle furnace at 550 °C (Wetzel & Likens 1991). By following these procedures, the fine (FIS) and coarse inorganic (CIS) sediment fractions in addition to fine (FOM) and coarse organic matter (COM) were determined.
To determine the Al, As, Cd, Cr, Cu, Ni, Pb and Zn concentrations in the sediment samples, the sample material was dried (40 °C) and then digested in triplicate with HNO 3 , HCl and H 2 O 2 on a digitally controlled heater plate (Quimis, Q313M) according to method 3050-B (US-EPA 1996). Readings were performed in an ICP-OES (Agilent Technologies 700 Series). Quantification limits were determined according to the equation: LQ = 10s/α.

Water
Analyses of the pH, dissolved oxygen (DO) and electric conductance (EC) in the bottom water were performed in situ by using a Horiba -U50 multiprobe. The water transparency was measured with a Secchi disk (SD). Samples were   (US-EPA 1992). Readings were performed in an ICP-OES, and the limit of quantification was calculated. Fortified solutions were prepared as checkpoints at concentrations of 0.20, 0.50 and 5.00 mg/L.

Toxicity
In this work, we considered total concentration of metals in the environment as a practical decision. We do not disagree that metal redox status, complexation or adsorption may interfere in toxicity (Hughes 2002, Flynn et al. 2014. However, such environmental conditions are not stable along the life time of BMI (Miao et al. 2006).
Considering that sediments where benthic macroinvertebrates live works as a pollutant sink, we believe that an approach considering total concentrations must give a more realistic situation of the ecosystem and the organisms experience, instead an instantaneous picture of a short-time situation. Moreover, many of these organisms ingest suspended matter and sediments (Cummins & Klug 1979, Merrit & Cummins 1996, Heino 2008. As consequence, BMI are subjected to toxicity by many pathways. Faunal composition must reflect the sum of distinct situations along time. The resulting metal concentrations were compared with toxicity reference values for ecological risk assessment (TRV) as calculated on the basis of the lowest effect level, LEL (Persaud et al. 1993) or effects range-low, ERL (Long et al. 1995) approaches depending on the analyzed metal according to US-EPA (1999) recommendations. The TRV considered here are based on toxicity assays that considered the total concentrations of metals. The ERL and LEL are thresholds at which toxic effects are expected to occur occasionally. At ERL, the lower 10 th percentile of observable effects in organisms will define the threshold. The LEL is the 5 th percentile of the screening level concentration (SLC). According to Persaud et al. 1993, "the SLC is an estimate of the highest concentration of a   contaminant that can be tolerated by a specific proportion of benthic species". The toxic units (TU) for each metal were then calculated with the equation TU = MC/TRV in accordance with Sprague (1970) and Castro-Catalá et al. (2015), with modifications.
Where MC is the concentration of the metal in the sediment or the water sample, and the TRV is the toxicity reference value for that metal. TU values that are higher than this unit indicate toxicity, and the more the TU exceeds this unit, the more toxic the sediment or water under analysis. The TRVs are shown in Table III.
The TRVs are based on US-EPA standards (1999). The TRV for metals in the water were determined by considering the dissolved fraction, with exceptions for Al and As, for which the total concentrations were considered.
TRV values considered here is just used as a reference and selection criteria for models building. They must not be considered as strictly threshold values in our work where total concentrations are measured instead of the dissolved ones. However, as the effects of dissolved concentrations and toxic units are dependent variables of total concentrations, the establishment of conceptual models must not be affected by such imprecision.

Statistical analysis
Pearson correlations (r) were calculated for the metal concentrations, and only metals that were not correlated could be entered into the models (r < 0.70). When a correlation was observed, the metal with the highest toxicity among all of the samples was chosen to use in the model to avoid autocorrelation.
Once the literature showed that the biological variables used here can respond to eutrophication, t-tests were performed between the dataset of biological variables against the total phosphorus (TP), chlorophyll-a (CL) and trophic state index (TSI) to check if alterations in these variables could influence the biological responses.
A group of environmental variables that were not directly related to metal toxicity are included, namely the pH, DO, EC, SD, FOM, COM, FIS and CIS. These measurements reflect the sediment composition, water physical-chemical conditions and trophic state in a broad sense. A Principal Component Analysis was performed by using these variables that were previously standardized as z-scores. This analysis was performed with the aim of including the environmental gradient, which is known here as "neutral", in the models with metal toxicities. To perform this analysis, the eigenvalues from principal component 1 (PC1) were considered. The scores of PC1 from this PCA is called here as ENV component. Biological data were resumed into five components representing the five biological levels on Table I plus abundance. Transformation procedures were like those described for ENV variables. As a result, models considered variations in biological levels (Species, Subfamily, Family, Tolerance, Function and Abundance) as a response to ENV and to metal toxicities.
The z-score standardization avoids the attribution of different weights based on the range of variation for each variable. With such procedure, the recorded value for each case (considered separately by variable) is transformed in terms of deviance of the mean. Considering this, z-score transformations remove the effect of greater abundance of some dominant species and permits the comparison in terms of variability and distribution.
Conceptual models were obtained by using a linear mixed effects approach. Linear mixed effects models include both random and fixed predictors, and they are recommended when working with block or nested designs. In the present work, each reservoir was considered as a block that could exert some random effects on a dependent variable (biological data), and the intrinsic variations that were promoted by the reservoir effects were therefore not neglected in the models (Zuur et al. 2009, Logan 2010. The models were selected by using information theory with Akaike Information Criteria (AIC). From a broad view, AIC measures the loss of information that occurs when choosing any given model in relation to the full reality. Models with smaller AIC values are then indicative of a higher fidelity to the full reality and are initially considered as the "best hypothesis". Important measures of the strength of a competing model are derived from the ΔAIC (the difference between the "best hypothesis" and any other competing model). One is the probability of a model (w, 0-1), which measures the strength of the evidence. Another is the evidence ratio (w/w 0 ), which provides information about "how many times" one model is better than another one. A comparison of the best model (w, the selected model) against a reduced model (w 0 , the model without any predictor variables but with only the intercept constant) will demonstrate how strong the selected model is. For more details about statistics based on information theory and criteria for model selection based on the I-T approach, see Burnham et al. (2011).
In the present work, we selected the six best models to explain variations in BMI indices and metrics grouped by levels by considering taxonomic measures with coarse and fine resolution in addition to non-taxonomic measures (functional, tolerance-based and abundance measures) based on sediment and bottom water characterization with an emphasis on the role of metal pollution.
For each biological variable listed in Table V, a full model containing the neutral gradient and all metals with a TU > 1 without correlations between them was generated. From this full model, the other 31 were generated and compared, including the reduced one. A total of 192 models were generated.
The selection of the best model for each biological level was based on the AIC and the importance (i v ) of the variables (i.e., the proportion at which one predictor variable was significant for explaining the descriptor variable by considering the 32 models that were generated) by using the following criteria: 1) Only variables with (i v ) ≥ 0.50 were considered.   2) The best model is preferably the one that includes all variables with (i v ) ≥ 0.50.
3) The evidence ratio that compares the best model according to criteria (1) and (2) with the reduced model must be at least 2.0 to be considered "significant". This criterion indicates that the selected model is at least twice as strong as a random distribution (the reduced model).
The mean (i v ) values were calculated for each predictor variable by considering the selected models as a general measure of the importance of each environmental component in the BMI assemblage.

BMI assemblage
In checking every sample, 36 species belonging to eleven subfamilies and seven families were identified. From this identification, only three species, two subfamilies and two families were collected from the three reservoirs. Chironomidae and Naididae were the most representative taxonomic families, and Chironominae, Tanypodinae (Chironomidae) and Tubificinae (Naididae) were the subfamilies that could be selected according to the occurrence criteria. The chironomids Polypedilum sp. and Tanytarsus sp and the tubificid L. hoffmeisteri were the most representative species.
There were up to five taxonomic subfamilies per sample. Only one subfamily was recorded in JC2, CA1 and CA2. The quantitative results were the same when comparing the number of subfamilies with the number of families per sample. Despite this finding, the composition was not the same (i.e., families comprised more than one subfamily, but in general, only one was distributed around all three reservoirs according to our selection criteria).
Polypedilum sp. was predominant in Cachoeira samples, and Tanytarsus spp. prevailed in Jaguari-Jacareí. L. hoffmeisteri has intermediary relative abundance in the Jaguari-Jacareí and Paiva Castro reservoirs. Lower density values were recorded in samples from Cachoeira (96 ind/m 2 ± 112). Paiva Castro (833 ± 362) and Jaguari-Jacareí (935 ± 398) were similar in relation to the density of organisms. A higher number of species occurred in Paiva Castro (11-15) than in Jaguari-Jacareí (5-10), and the lowest values were recorded in the Cachoeira samples (1-7). The species dominance followed a gradient of PC < JG-JC < CA. Taxonomic composition and abundance are showed in Table IV. The functional feeding groups known as collector-filters (mean = 30% ± 35), predators (27% ± 35) and collector-gatherers (26% ± 29) were predominant, and they had similar relative abundances per sampling point. However, their distribution was distinct. The predator guild was predominant at JG-JC, the collector-filterers dominated CA and the collector-gatherers were better represented in PC. A higher number of functional groups was observed in the PC samples (Table V) (2012). According to the index, all sampling points were under "regular" conditions.

Environmental characterization and toxicity
The bottom water from the reservoirs that was analyzed in this work presented a neutral pH (6.88±0.52), satisfactory oxygenation (8.99 ± 1.54 mg/L) and a warm temperature (17.99 ± 0.59 °C). The Secchi disc readings were approximately 2.14 ± 0.63 m, and moderate electric conductivity was recorded (35.55 ± 5.34 µS/cm). The sediments were made up of more than 10% organic matter at JC1, CA1, PC1 and PC2. The organic matter at Paiva Castro was primarily coarse relative to that at JG-JC and CA. JC2, JC1 and CA2 were the sampling points with higher coarse inorganic sediments (62.26 ± 0.59%). The Paiva Castro inorganic sediment was primarily formed by fine particles (82.55 ± 8.31%).
Principal component 1 (ENV) explained 62.63% (PCA1 broken stick = 34.73%) and could be used to distinguish among the three reservoirs. An increase in the ENV occurred when the COM and FIS decreased and when the CIS and pH increased (Figure 3).
The metals Al, As, Cd, Cu and Cr had a TU ≥ 1 in one or more sediment samples. In addition, Al, Cr and Cu had a TU ≥ 1 in the bottom water. The As, Cd, Ni and Pb concentrations in the water could not be determined once they fell below the limit of quantification (Table VI).
Correlations were observed between the metal toxicity in the sediments and some other environmental variables, with As and Al in relation to FOM (r = 0.80 and 0.88), Cu and Cd with COM (0.83 and 0.91), Cu_s with FIS and CIS (0.66 and -0.71) and Cr with FIS (-0.69) and CIS (0.65).
A correlation between metal toxicities in sediments was observed between As-Cu (0.72), As-Cd (0.75) and As-Al (0.84). As a consequence, data about As toxicity were not considered in the models (it had a lower toxicity than Cu). Because Cd was toxic in only one sample, it was not considered either. With regards to the water analysis, only Al had a TU ≥ 1 in more than one sample, and it was not correlated with any other variable.

Biological and environmental interactions and the linear mixed effects models
By comparing trophic state data with biological variables, it could be concluded that the trophic conditions did not influence the BMI metrics and indices in our research (Table VII). Some biological variables were correlated. In such cases, the simplest one was maintained, and the correlated variable removed to avoid double weight in the biotic PCA. Dominance (correlated with Polypedilum sp. and species richness); DF (correlated with CHD, NAI and FAM); SF (correlated with TNP); Hilsenhoff index (correlated with TOL and TNT: CHD); and D.FUNC (correlated with FUNC) were previously removed and did not entered to the biotic PCA analyses.
The PCA-based variables can be broadly described as follows (Figure 4): -Species/Genus level: low values indicate the predominance of L. hoffmeisteri while high values point to higher abundances of Tanytarsus spp.
-Subfamily level: low values were weakly associated with higher number of subfamilies while higher values were indicative of higher abundances of Chironominae (weakly) and Tubificinae (strong relation).
-Family level: low values were weakly associated with high number of families while high values were strongly associated with higher abundances of Chironomidae and Naididae.
Function level: higher values were indicative of higher abundance of collector-gatherers.
-Tolerance level: higher values were weakly associated with higher abundance of organisms considered as "tolerant" and strongly associated with higher TNT: CHD rates.
-Abundance: PC1 values were related to lower density (log 2 ) values.
In considering all of the selected models, the most frequent abiotic variables were the toxicity of aluminum in the water (Al_w) and the toxicity if chromium (Cr_s) in sediments (importance ≥ 0.50). The neutral environmental gradient and other metal toxicities could not be considered as good predictors of the biological data. The biological variable "tolerance" did not respond to any group of abiotic factors considered in this research. The strongest model was that described by "subfamily" variable (2417.53 times stronger then the reduced model).

Environmental characterization and toxicity
The sediment grain size, organic matter content and pH were relevant to the BMI assemblages, and they presented a well-established gradient along the reservoirs. Previous studies indicated the relevance of some of these environmental   variables to the composition and distribution of BMI (Beghelli et al. 2012, Jones et al. 2012. In addition, changes in these environmental variables tend to alter the distribution of the metals in the environment, as well as their bioavailability , He et al. 2012, Rangel et al. 2011). The correlations obtained in the present study corroborate this assumption, indicating the interrelations between the sediment composition (the predominant grain size and coarse organic matter proportions) and metal concentrations in sediments. The proportion of COM, CIS and FIS in sediments may be related to the integrity of the riparian forest, which has an important role in the sediment composition from falling debris and acts as a barrier to the movement of particles from land to water (Naiman & Décamps 1997).
In considering the sediment toxicity, the following two primary groups could be distinguished: a low toxicity group formed by PC1, PC3, CA2 and CA3 and a high toxicity group that includes Jaguari-Jacareí sampling points, PC2 and CA1. Despite the fact that PC2 could be grouped with JC-JG and CA1 in relation to general toxicity, it must be qualitatively distinguished. The toxicity of PC2 is related to Cu_s and Cd_s, and the primary contaminants in JC-JG and CA1 are Al and Cr.
The water toxicity was primarily determined by the Al, which was present in higher concentrations than the established TRV in all of the samples. Despite a well-documented effect of acidity on Al toxicity related to its low solubility when in neutral conditions, many works with benthic macroinvertebrates demonstrated that even in such cases aluminium is still bioavailable to BMI (Michahilova et al. 2003, Shuhaimi-Othman et al. 2011. The toxicity from aluminum in the water varied from 1.6 in CA2 (0.139 mg/L) to 177.4 TU in JC1 (15.433 mg/L), indicating values that were many times greater than the LEL or ERL for test organisms as established by the US-EPA 1999.
The highest values were concentrated upstream where there is a lower proportion of natural The specific situation of the PC2 sampling point must be related to the water flow, macrophytes and riparian forest and also to human impacts. This sampling point is in a backwater area where the presence of macrophytes and a denser riparian forest could be observed. This situation favors the reduction of water flow and sedimentation, promoting metal and organic matter accumulation (Martins et al. 2015). The positive correlation of Cu and Cd with COM and Cu with FIS are indicative of the tendencies of these metals to adsorb to these sediment fractions, which were especially high in Paiva Castro (Demirak et al. 2012, Martins et al. 2015. Despite having the highest proportion of natural vegetation, this area also has the highest proportion of urbanized areas, and pollution by multiple metals is common in such cases (Zhang & Liu 2014).

Biological responses
It is necessary to mention here that seasonality effects were not considered in the present work. Instead, we optioned to replicate sampling design during the same season to avoid an extra source of variability. Effects of seasonality in benthic macroinvertebrates composition, are usually associated with changes in organic matter, dissolved oxygen or in trophic state conditions. These variables were monitored here (Linares et al. 2013, Beghelli et al. 2014.
Given the complex differences the reason to avoid seasonality effects is similar to those that drives laboratorial assays: we lose some environmental information to get results that are more precise by controlling environmental variability.
Considering that the aim of this research was to get effects of toxicity by metals on BMI by generating conceptual models of type doseresponse, the results and conclusions here are not expected to be significantly affected by seasonality.
T-tests on the trophic state showed that the trophic conditions were not affecting the biological variables studied here. As mentioned, many works have already presented some of the metrics and indices considered in this research (the TOL, HIS, the proportion of Limnodrilus or tubificids, species richness, and dominance) to eutrophication or a "general" environmental condition (Beghelli et al. 2014, Qin et al. 2014). This result is likely more a consequence of low variations in the trophic state among sampling points than the true absence of a relation.
Furthermore, in a broad sense, eutrophication is made up of other changes that are more directly related to BMI assemblages than the surface TP, TN or chlorophyll-a concentrations because increasing the organic matter in sediments depletes the dissolved oxygen on the bottom, in addition to pH alterations and the death of macrophytes (Takamura et al. 2009, Nelson and Steinman 2013, Qin et al. 2013). The first three factors were considered in the neutral environmental gradient (ENV). The BMI assemblages responded to metal toxicity in water from Al and in sediments from Al, Cr and Cu. Moreover, the ENV variable consisting of the COM, FIS, CIS and pH was relevant and helped to explain seven of the 25 biological variables.
Novelties and differences between the present work and other similar studies should be highlighted. Zhang & Liu 2014 performed a comprehensive study on BMI assemblages and periphyton responses to the environment. Other authors used an ecological risk approach to study human impacts on three distinct habitat types. These authors related the correlation of BMI metrics to the ecological risk index for metals in distinct areas. In comparison with the present research, the previous study has the advantage of being more comprehensive by including many stressing conditions in addition to two distinct biotic assemblages. However, our research presents some advances in the sense that the present study is more specific. First, the work of Zhang & Liu (2014) presented many sources of variation such as seasonality and organic pollution, and furthermore, their pattern of environmental conditions (toxicity from distinct metallic elements presented the same pattern) and their statistical approach would not permit the differentiation of specific metal effects against other environmental variables.
The present work highlighted more specific interactions by considering the toxicity of distinct elements in sediments and water in addition to the random effects of reservoirs and other relevant environmental characteristics in the ENV component. Furthermore, the toxicity caused by individual metals in the Cantareira Complex did not present similar patterns to one another, and thus, the toxicity from Al_w, Al_s, Cu_s and Cr_s could be clearly differentiated. The trophic state did not play an important role in our biological data.
The biological variables that were analyzed in this work were based on proportions and not on direct abundance, with the exception of the density and taxa numbers. Under these considerations, it can be assumed that any positive relation between a specific taxonomic or functional group and metal toxicity is not the consequence of a benefit that was directly caused by toxicity on the biological variable. Instead, it is much more plausible that it is a consequence of a proportional reduction in other more sensitive groups ((Mackintosh et al. 2015, Milner et al. 2016.
Only Al_w and Cr_s could be considered as good predictor variables to BMI assemblage composition and distribution. Variations among taxonomic and functional indexes and metrics could be attributed to variations in sensitivity to specific metals considering characteristics such as escape strategies, physiological mechanisms of toxicity avoidance, the degree of exposition to sediments or water and feeding habits (Méndez-Fernández et al. 2013, Clements et al. 2015, He et al. 2015. Environments with high Cr_s and low Al_w concentrations could be associated to predominance of Chironomidae and Naididae when considering the family level and Chironominae and Tubificinae (subfamily level). When species/genus level is achieved, we predict high relative abundances of Tanytarsus sp. (Chironominae) in chromium rich sediments (low Al_w) and high relative abundances of L. hoffmeisteri (Tubificinae) in environments with high Al_w values (low Cr_s). It is important to notice that only the finest taxonomic resolution level could detect this difference among chironomids and oligochaetes.
As we expected, tolerance indices were not good predictors to metal pollution in our study. The main explanation for that must be the relation of such classification of "tolerance" based on organic pollution, eutrophication or deforestation. In environments where toxicity is not associated to organic pollution or eutrophication, these indices may vary in a different way. Similar results were observed by López-Dovál et al. 2012, who tested the response of an index on the basis of general pollution tolerance that was similar to that of HIS against the TU for chemicals, especially metals, in four Spanish basins. The authors did not find significant correlations between the variables.
Despite to has some power of prediction, the functional feeding guilds were the weaker model obtained been just 2.71 times more probable than the reduced. Composition in terms of food intake are related essentially to the entrance of organic matter from the surroundings and to the water flow. The relation of such kind of biological traits to metal pollution must be consequence of metal binding to organic matter and sediments instead toxicity.
Our model suggests the occurrence of an assemblage dominated by collector-gatherers as the concentration of chromium in sediments increases. Considering that toxicity by Cr_s was higher in reservoir with less natural vegetation some care must be used to interpret this model. Environments with higher proportion of riparian forest are expected to shelter higher diversity of feeding groups than environments with other land uses (Cummins et al. 2005, Compin & Céréghino 2007.
In fact, ecotoxicological assays have demonstrated toxic effects of Al_w on BMI organisms even at neutral pH. Results reported by Shuhaimi-Othman et al. 2011 are representative of the two predominant families and subfamilies recorded in this work.
Concentrations higher than LC 50 reported for C. javanus (Shuhaimi-Othman et al. 2011) were recorded in four cases (JC1, JC2, CA1 and CA3) while an Al_w concentration higher than the LC 50 value reported for N. elegans were recorded in sample from JC1.
Density decreased as toxicity by aluminium in water grew. Considering the broad depressor effect of toxicity on biota and that Al_w was the toxicant with highest values in terms of TU, this result reinforces the idea that assemblage abundance measures can be suitable to assess the most severe impacts related to toxicity in a specific environment.
The taxonomic approach was clearly the group with the best predictors to metal pollution in Cantareira Complex. The models obtained here follow a gradient of strength: family < species/ genus < subfamily taxonomic level. Despite many debates about the ideal taxonomic resolution for bioassessment purposes, researchers usually agree that the choice of the best resolution must consider some compromise among scientific ideal, bioassessment purposes, financial and logistic realities (Jones 2008, Buss & Vitorino 2010, Frizzera & Alves 2012. Our results show that subfamily level of resolution was the most adequate to monitoring toxic effects of metals in Cantareira Complex with the lower loss of truly valid information to predict toxicity. It means that when taxonomic resolution was limited to family level, many information was lost (the model was 76.75 times weaker) in comparison to subfamily level. However, when comparing the subfamily level with the finest resolution, the subfamily could be considered as the most adequate and parsimonious choice.
In considering our results and those in the literature, we strongly recommend the development of specific tolerance indices for metal pollution in future studies on BMI assemblages. In this respect, the present work provided important data by demonstrating a set of indices and metrics that truly respond to metal toxicity as also compared distinct taxonomic resolutions. One advance made by this work was that the trophic state did not present significant effects, and other environmental components that were frequently relevant to BMI assemblages were analyzed in the models to assess the effects that were caused separately by metals.
The elaboration of a general toxicity index may present some problems related to the differential sensibilities of organisms to specific toxicants. Indices based on total tolerance scores may fail exactly because of these specific differences. For instance, L. hoffmeisteri is generally considered as "tolerant", but it is more sensitive to chromium contamination in sediments than Tanytarsus sp. which is usually considered as a "sensitive" organism.

CONCLUSIONS
The present work showed the effects of metal toxicity (potential risk) in water and sediments on BMI assemblages. The biological variables responded to the toxicity at different levels (coarse taxonomic resolution, fine taxonomic resolution, functional indices and density).
From our results, we concluded that: the variables from the taxonomic group were better bioindicators. Subfamily taxonomic resolution was the most adequate to predict metal toxicity. The functional feeding group responses to metal toxicities may be related to metal adsorption and complexation with organic matter and sediment grains as also to the food ingestion and mobility of organisms.
With our findings, we provided the basis to develop a new BMI index for metal pollution in oligotrophic reservoirs in tropical environments. Further studies may apply similar approaches in eutrophic reservoirs. The application of these