In situ toxicity and ecological risk assessment of agro-pesticide runoff in the Madre de Dios River in Costa Rica

The River Madre de Dios (RMD) and its lagoon is a biodiversity rich watershed formed by a system of streams, rivers, channels, and a coastal lagoon communicating with the Caribbean Sea. This basin sustains a large area of agricultural activity (mostly banana, rice, and pineapple) with intensive use of pesticides, continually detected in water samples. We investigated in situ the toxicological effects caused by pesticide runoff from agriculture and the relation of pesticide concentrations with different biological organization levels: early responses in fish biomarkers (sub-organismal), acute toxicity to Daphnia magna (organismal), and aquatic macroinvertebrate community structure. The evaluation was carried out between October 2011 and November 2012 at five sites along the RMD influenced by agricultural discharges and a reference site in a stream outside the RMD that receives less pesticides. Acute toxicity to D. magna was observed only once in a sample from the RMD (Caño Azul); the index of biomarker responses in fish exposed in situ was higher than controls at the same site and at the RMD-Freeman. However, only macroinvertebrates were statistically related to the presence of pesticides, combined with both physical-chemical parameters and habitat degradation. All three groups of variables determined the distribution of macroinvertebrate taxa through the study sites.


Introduction
Costa Rica's economy is largely based on agriculture though ecotourism is rapidly gaining ground due to the country's high biodiversity and pristine nature (ICT 2016). This creates a conflict of interest between the deleterious impact of agriculture and agrochemicals on the environment and the need to preserve its biodiversity and natural habitats.
The River Madre de Dios (RMD) watershed is located in the Caribbean region of Costa Rica. It includes a network of streams, creeks, rivers, and channels, which drain into a coastal lagoon connected to the Caribbean Sea. Globally, coastal ecosystems are characterized by their high biodiversity and constitute a very valuable habitat for marine, estuarine, and freshwater species, while also sustaining important economic activities (Newton et al. 2013;Pérez-Ruzafa and Marcos 2012). Besides their ecologic and economic relevance, estuarine and freshwater coastal environments and their ecosystems face the threat of pollution and other anthropogenic pressures (Hernández-Romero et al. 2004;Kuk-Dzul et al. 2012). In the case of the RMD, 52 % of the basin is dedicated to agriculture (mainly banana, pineapple, and rice) with an intensive use of pesticides that pose a risk of contamination to this important ecosystem (Bravo et al. 2013;Diepens et al. 2014).

Responsible editor: Ester Heath
The amount of pesticides used in banana in the country reaches up to 76 kg of active ingredients (a.i.) per hectare and year (kg a.i./ha/year) mainly as fungicides, nematicides-insecticides, and herbicides. In pineapple cultivation, 47 kg a.i./ha/year is used mainly in formulations of insecticides-nematicides, herbicides, and fungicides. In rice fields, 7 kg a.i./ha/year is used in the form of herbicides, fungicides, and insecticides . These crops are widely distributed in the Caribbean region, where studies have found high concentrations of pesticides downstream from agricultural activities in rivers in the area (Castillo et al. 2000Echeverría-Sáenz et al. 2012). Detection of such substances has been also associated with massive kills of aquatic fauna (crustaceans, reptiles, birds, and mostly fish) throughout the Caribbean, where pesticide residues have been detected and quantified in water and biota at the same moment of the event (CGR 2013; Polidoro and Morra 2016).
For more than a decade, chemical pollution has been observed in the Madre de Dios lagoon and its tributaries (Diepens et al. 2014;Mena et al. 2014a). Efforts to characterize toxic effects of water and sediment samples from the RMD have been carried out with laboratory toxicity tests using fish, zooplankton, and benthic invertebrates (Daphnia magna and Hyalella azteca) (Ålander 2006;Diepens et al. 2014). Also, more recently, Arias-Andrés et al. (2016) have estimated the toxic risks of individual and mixtures of pesticide residues found in the surface water of the RMD from 2009 to 2012, using the species sensitivity distribution, msPAF, and PERPEST model approaches.
In order to assess whether risks related to the presence of pesticides translate into effects in the ecosystem, an evaluation of effects at different levels of biological organization was performed in the RMD. A combination of methods that assess early warning signals and higher level ecological effects of pollution would give a more integrated overview of the environmental health in an ecosystem (Vasseur and Cossu-Leguille 2003;Peschke et al. 2014;Martinez-Haro et al. 2015). The implementation of in situ approaches, such as exposure of organisms in the field, allows for the observation of effects in a real scenario where an interaction of environmental parameters and pollutants occurs (Crane et al. 2007).
Given this context, this study aims to investigate the effects caused by pesticide runoff from agriculture to the RMD, at three levels of biological organization: (a) biochemical level (by measuring a set of biomarkers in native fish exposed in situ), (b) individual level (acute toxicity of surface water to D. magna), and (c) community level (monitoring macroinvertebrate community structure) in the rivers and canals of the RMD.

Materials and methods
The field study was carried in the months of October-November 2011 and February-March and October-November 2012. December through January and May through August were avoided due to very heavy rainfalls, which could increase the probability of losing experimental equipment or not being able to access the sampling points due to frequent river overflowing. Four sites along the RMD basin were sampled: site 2 (CA-S) on Caño Azul stream, a tributary of the RMD that drains banana and pineapple plantations; site 4 (CPama-J), located at the Pama Canal which receives water from a banana plantation; site 5 (RMD-F), in the Madre de Dios River before joining the Pama Canal; and site 6 (URMD-CPama), after the junction of the Pama Canal and RMD. Site 13 (EARTH) was selected as a reference site in a stream outside the RMD, close to a banana plantation which is managed by an agronomy university (EARTH) and is less exposed to pesticides and fertilizers than the RMD, due to better agricultural practices (www.earth.ac.cr/en/aboutearth/earth-products/banano/). Also, it is slightly upstream of the banana plantation and its drainage canals (Fig. 1, see Table 1 for a detailed description of the sampling sites). Surface water samples were collected twice each sampling occasion with a 48-h interval between them (first, when fish cages were deployed, and second, when they were retrieved). Water samples were collected by inserting pre-washed 2-L brown glass bottles into the water. Temperature (°C), pH, dissolved oxygen (mg/L), and conductivity (μS/cm) were measured in situ with a portable multi-probe (Hach HQ40 d). The collected water samples were transported in refrigerated coolers to the laboratory (LAREP, UNA) and stored at 4°C in the lab for less than 24 h before analyses of nutrients, suspended solids, and pesticide residues.
Determination of pesticide residues, nutrients, and total suspended solids Pesticide residues in water samples were analyzed using gas chromatography with mass spectrometry (GC-MS) for nonpolar pesticides and high-performance liquid chromatography (HPLC) with diode array detection for polar pesticides (LC-PDA). Nutrients (NO 3 , total and soluble phosphorus) and total suspended solids (TSSs) were determined according to the Standard Methods for the Examination of Water and Wastewater (APHA 2005). Further details on chemical analyses and a thorough listing of pesticides used in the studied area and analyzed in the laboratory are included in Castillo et al. (2016).

Conversion of pesticide concentrations to toxic units
Measured concentrations of pesticides were transformed into toxic units (TUs) for arthropods, using data on effects collected from the ECOTOX-EPA database (US EPA 2015), the Pesticide Action Network database (PAN 2015a), the E-toxBase database (De Zwart 2002), as well as data previously reported by studies carried out at IRET with native species. Data were filtered, and only freshwater laboratory experiments that evaluated effects on immobility, mortality, and population endpoints in short-term experiments (1-4 days) resulting in an EC 50 or LC 50 were included. TU values for each pesticide (with available information) were based on the median effect concentration for arthropods and calculated according to the following formula:  (1.05 × 108 cells/L), YCT (yeast, cereal leaves, TetraMin; 0.5 mL/L), vitamin B 12 (2 μg/L), and selenium (2 μg/L) (US EPA 1982, 1989WDNR 1991;Mount and Mount 1992). For the assay, ten neonates younger than 24 h old were exposed in triplicate to 30 mL of the samples and kept at 20 ± 2°C. Hard reconstituted water (HRW; Dutka 1989) was used as a control. Observations of organisms' mobility were made at 24 and 48 h of exposure.

In situ test with caged fish and biomarker analyses
Immature juveniles of Parachromis dovii (mean standard length 2.98 ± 0.29 cm) were acquired from the Aquaculture Laboratory at the Universidad Nacional, Heredia, Costa Rica. The fish were transported to the field in aerated containers filled with clean culture water and then deployed in cages at each site. Rectangular cages were built using a PVC pipe frame (25 cm × 30 cm × 15 cm) covered with a 0.4-cm mesh size net. Four fish were added per cage, using two cages per site. In each test, a control (non-exposed) group of fish (n = 8) was kept in aerated containers with clean culture water on land until the end of exposure. After a 48-h period, caged and control fish were dissected and samples of their liver, brain, and dorsal muscle were collected, frozen in liquid N 2 , and preserved at −80°C, until biomarker analyses. This 48-h period was adjusted after attempts of carrying a 96-h exposure and losing cages and fish due to abrupt changes in the levels and currents at the sampling sites, which are frequent in tropical rivers. Besides, there are also reports of significant responses of biomarkers in fish after exposures of 48 h to pesticides in laboratory assays or in caging experiments (Sayeed et al. 2003;Mena et al. 2014b).
Muscle and brain cholinesterase (ChE) activity as well as glutathione S-transferase (GST) activity, lipid peroxidation (LPO), and catalase (CAT) activity in the liver were evaluated. Samples were homogenized in an appropriate buffer and processed as described in Mena et al. (2014b) for biomarker determinations. Briefly, protein content in sample homogenates was determined by the method of Bradford (1976) using γ-globulin as a standard. ChE activity was measured using the method of Ellman et al. (1961), using 1 mM acetylthiocholine and 0.1 mM 5,5′-dithiobis-2-dinitrobenzoic acid (DTNB) as a substrate and conjugate, respectively; reaction was measured at 415 nm during 15 min and expressed as nanomoles per minute per milligram protein. GST activity was determined as described by Habig et al. (1974), exposing samples to 1 mM CDNB and 1 mM GSH and monitoring the reaction at 340 nm during 3 min; activity was reported as nanomoles per minute per milligram protein. Lipid peroxidation was measured by the thiobarbituric acid reactive species (TBARS) assay (Oakes and Van Der Kraak 2003) and expressed as nanomoles of TBARS per milligram of protein.
CAT activity was measured according to Aebi (1974) by the decrease in absorbance at 240 nm during 20 s due to H 2 O 2 consumption and expressed as micromoles per minute per milligram protein. Biomarker data were transformed to generate an integrated biomarker response (IBR) according to Devin et al. (2014). As suggested by these authors, circular permutations in the order of biomarkers were created during the calculation to generate a set of IBR values for every sample which allowed the comparison between sites and control.

Macroinvertebrate community structure
Four artificial substrates (clay bricks 20 × 5 × 5 cm, placed inside plastic mesh bags, with 0.3 cm mesh size) were placed at each sampling site and left to be colonized by benthic invertebrates for 1 month. Artificial substrates were tied to branches of trees in the river bank and submerged no more than 50 cm into the water. The weight of the brick itself submerged the artificial substrates, and the depth was controlled simply by adjusting the length of the cord. Samples were recovered after 1 month, filtered through a 250-μm sieve, and preserved in plastic containers with 96 % ethanol. Organisms were sorted in the laboratory under a stereoscope and then identified to family and genus levels (when possible) using relevant taxonomical keys (Springer et al. 2010;Merritt et al. 2008;Pennak 1978).

Statistical analyses
In order to assess the significance of the correlations between the community composition of macroinvertebrates and the three groups of explanatory variables (pesticides, habitat parameters, and physical-chemical variables), Monte Carlo permutation tests were performed using redundancy analysis (RDA). The significance of each explanatory variable was assessed separately. Before the analysis, the abundance counts and pesticide concentrations were Ln(2x + 1) and Ln(200x + 1) transformed, respectively, according to Van den Brink et al. (2000). Missing values were replaced by the median value of the measured values. Permutations were performed by permuting within sampling dates only. A separate RDA was performed with all significant explanatory variables to construct a biplot showing the correlations between the explanatory variables, sites, and species.
After this, variance partitioning of the macroinvertebrate data was performed between the three groups of explanatory variables, i.e., pesticides, habitat parameters, and physicalchemical variables. To limit the number of explanatory variables, only those evaluated as significant by the earlier performed Monte Carlo permutation tests were included. This analysis provides an overview of the unique part of the variance that is captured by the three groups and which amount is shared between them.
For the biomarker dataset, similar Monte Carlo permutation tests were performed using only the pesticides and physical-chemical parameters as explanatory variables. The biomarkers were standardized and centered within the analysis, making them mathematically equally important irrespective of their measuring scale.
Differences in IBR values were assessed with a Kruskal-Wallis test coupled to a pairwise comparison using IBM ® SPSS ® Statistics 22 software. All other analyses were performed using Canoco 5 (Ter Braak and Šmilauer 2012). For more background information on the use of ordination for the analysis of ecotoxicological datasets, see Van den Brink et al. (2003).

Results and discussion
Site characterization and physical-chemical properties The assessed sites (Table 1) differ in certain characteristics that could impact their exposure to pesticides, such as relative proximity to cultivation fields. For instance, sites 2 and 4 in the RMD, as well as reference site 13, are closer to crops compared to the other two sites in the RMD (sites 5 and 6). All sites receive runoff from banana plantations, but sites 2 and 4 are also influenced by rice and pineapple plantations that use other types of pesticides which can also be transported downstream to sites 5 and 6 ( Fig. 1).
The riparian vegetation is far more conserved at site 6, where the water course is practically in its natural state, while upstream, sites 4 and 5 present some degree of deforestation, with pastures along one side of the river bank. The Caño Azul stream (site 2) and the reference site EARTH (site 13) only have a few meters of vegetation cover, mostly composed of patchily distributed trees (site 2) or shrubs (site 13). These sites are also shallow (<1 m depth) stream-like water bodies, with moderate current velocity, while the other three sites are deep (>3 m), with low-velocity canal/lowland river types.
Physical-chemical parameters (Table 2) were within ranges appropriate for sustaining aquatic biota in four out of the five sites. CPama-J (site 4) was the exception with dissolved oxygen (DO) values just at the limit of 2 mg/L required for many aquatic organisms (Chapman 1996;Levin et al. 2009). Further, this site had the higher values of soluble phosphorus, both parameters related to events of eutrophication which is known to affect biodiversity, when oxygen depletion levels reach such low values (Chapman 1996;Correl 1998). Another factor which may have decreased oxygen concentration at this site was diminished water movement, since this canal was almost stagnant (0 to 0.6 m/s). Higher conductivities measured at sites 4 and 6 are also an indication of erosion that might be related to runoff. On the other hand, our reference site (site 13) showed physical-chemical parameters favorable for aquatic life during the whole sampling period. As site 13 is outside the RMD basin area, it slightly differs in elevation, climate, and ecological characteristics compared to the other sites.

Pesticide residues and TUs
Throughout the study period, 27 substances (13 fungicides, 7 herbicides, and 7 insecticides) were detected in surface water samples, including terbufos-sulfone (terbufos metabolite) (Fig. 2). The site with the highest average number of pesticide detections per sampling occasion during the study was site 2 (10.2 a.i.), followed by site 5 (9.8), site 6 (8.3), and site 4 (6.8). At the reference site (site 13), pesticide residues were detected only in three sampling occasions (chlorpyrifos in March 2012, terbutryn in October 2012, and epoxiconazole in November 2012). Maximum amounts of pesticide residues of each substance detected during the whole study are presented in Fig. 2. Transformation of the pesticide residues into TUs for arthropods showed that all sites in the Madre de Dios watershed carry a constant load of toxic pesticides, the highest levels being found in sites 2, 5, and 6; with the highest in site 2 during the October 2012 sampling (Table 3).
The majority of the pesticides detected were fungicides which are mainly applied in banana plantations. High concentrations of herbicides like diuron and insecticides such as diazinon in surface water could be related to the high use of these compounds in crops like rice and pineapple. Twelve (43 %) of the detected a.i. are classified as highly hazardous pesticides according to PAN (2015b). Arias-Andrés et al. (2016) performed a low-tier assessment based on the evaluation of over 200 water samples collected at the RMD watershed between 2009 and 2012. The mentioned study, which includes native species in the evaluation, calculated SSD for each pesticide and concluded that the higher risk to freshwater biota was mainly posed by herbicides (diuron, ametryn, bromacil, hexazinone, oxifluorfen) and insecticides (carbofuran, diazinon, fenamiphos, ethoprophos, and terbufos) constantly present in the water.
Further analyses by Rämö et al. (2016) indicate a high probability of direct toxic effects by a concentration of certain herbicides (ametryn, diuron, and bromacil) and insecticides (carbaryl, diazinon, ethoprophos, and terbufos), and higher risks for primary producers and arthropods compared to fish. Therefore, given that effects on arthropods have been studied in the present investigation, effects on primary producers should be assessed in following studies.
It should be considered that due to analytical limitations, some relevant pesticides were not analyzed in water samples. An example is mancozeb, and this dithiocarbamate fungicide is the most applied a.i. in banana plantations in the study area (26.1 kg a.i./ha/year) (Bravo et al. 2013). This compound is applied by airplane, and high levels of its metabolite ethylene thiourea (ETU) have been detected in human populations living near banana plantations in the area (Mora et al. 2014).

Laboratory toxicity tests
Acute toxicity (60 % mortality in an undiluted sample) to D. magna was observed only in organisms exposed to a sample from site 2 in one sampling occasion (November 2011). Mortality in HRW (negative control) was <10 % in all tests. Data is average (X) and standard deviation (SD) of 12 measurements done per site  Epoxiconazole was present in the sample (0.05 μg/L), but TUs were not calculated due to a lack of toxicity data that met our selection criteria Even though there were samplings with higher pesticide content and TU estimation during the study, this specific site and date had more TSS (23 mg/L) than the majority of samplings (median = 10.2 mg/L). According to Capper (2006), Levine et al. (2005), Kirk and Gilbert (1990), and McCabe and O'Brien (1983), suspended solids can be toxic to cladocerans, since they can ingest particles of sediment that subsequently become lodged in the gut tract, causing not only mortality but also decreased ingestion rates and assimilation efficiency, especially in juveniles. This can partly explain our results; however, there were other two samples on the same date at sites 5 and 6 with similar TSS concentrations which did not lead to Daphnia mortality. Another possible explanation would be that the suspended solids had affected the bioavailability of pollutants. Yang et al. (2006) found that the source and properties of SS could increase the adsorption and reduce the bioavailability of pyrethroids. However, we did not analyze the properties of SS in this study.

In situ fish exposure and biomarkers
Regarding the biomarkers measured in fish, ChE activity inhibition is probably one of the better characterized biochemical biomarkers used in ecotoxicology, indicating the presence and effects of organophosphate and carbamate pesticides on the nervous system (Bocquene and Galgani 1998). GST, LPO, and CAT were measured as biotransformation-and oxidative stress-related markers. These are more unspecific biochemical responses that can be altered by different pollutants, including pesticides (Hamed and Farid 2003;Valavanidis et al. 2006). After testing the association of the biomarkers with the physical-chemical parameters and presence of pesticides in the samples, only one variable explained a significant part of the variation in biomarker measurements (total phosphorus, p value 0.019). This was not expected, and we did not find references with information that could support this finding. All other parameters obtained a p value larger >0.1. Integration of biomarker data into IBR showed the highest magnitudes in biomarker responses in October 2011 (site 6), March 2012 (site 5), and February 2012 (site 5) while IBR values significantly higher than controls were more frequent in site 2 (Table 4). However, as no significant correlation between individual biomarkers or IBR and presence of pesticides could be established, their signals should be interpreted carefully. Considering this, efforts should be undertaken in order to verify cause-effect relationships between biomarkers and relevant pesticides and to identify other sources of variability. This tool coupled with other toxicological and ecological evaluations should contribute to the establishment of links between the presence of pollutants and effects at different levels of biological organization (Sanchez and Porcher 2009;Barata et al. 2007;Colin et al. 2016).
Another possible interference factor is related to pesticides (such as the fungicide mancozeb) which were not analyzed in water but have a high potential of being present in the waters of the RMD due to its frequency of application and the fact that it is aerially applied, therefore increasing the drift potential outside target areas. Although mancozeb's acute toxicity to fish is relatively low (US EPA 2015) and therefore fish are expected to survive in water with the presence of the substance, neurotoxicity (Jarrad et al. 2004) as well oxidative stress (Kubrak et al. 2012) biomarkers have been reported to respond to high concentrations of mancozeb.

Macroinvertebrate community structure
The analysis of macroinvertebrate community (MC) is used in this study as a higher level of biological organization endpoint, which, hand by hand with sub-individual-level (biomarkers) and individual-level (Daphnia toxicity test) endpoints, gives a wider array of information on the effects of environmental concentrations of pesticides on this ecosystem's aquatic biota.
More than 50,000 individuals from 13 classes, 21 orders, and 47 families were collected and identified (see the most abundant taxa in Table 5). The macroinvertebrate community structure in the reference site (site 13) outside of the RMD was clearly different from all sites within the RMD. This site had the highest diversity, the most sensitive species (representatives of Elmidae, Ephemeroptera, Trichoptera, Gomphidae), and the lowest pesticide concentrations. Amongst the other four sites, the macroinvertebrate community had the lowest species richness at site 4, which had the lowest diversity and the lowest DO of all (∼2 mg/L). Sites 5 and 6 had very similar community compositions, while site 2 shared species with all the other sites. These polluted sites had lower richness and also higher abundances of tolerant taxa (to organic pollution) such as annelids, molluscs, and dipterans (MINAE-S 2007; Table 5). The individual variables explaining more than 20 % of the variations in the structure of the macroinvertebrate communities were habitat parameters such as type of water body and vegetation on the river bank. Furthermore, variables like pH, DO, conductivity, and soluble phosphorus, as well as pesticides (pyrimethanil, diuron, ethoprophos, epoxiconazole, azoxystrobin, ametryn, diazinon, and metalaxyl), also explained a significant percentage of the variation (Table 6; Fig. 3).
Habitat variables are well-known influencing factors for distribution of macroinvertebrates (Jáimez-Cuéllar et al. 2002;Lorion and Kennedy 2009) and follow theoretical notions such as the river continuum concept (Vannote et al. 1980). Furthermore, physical-chemical parameters that reflect water quality, especially those related to organic pollution (DO, nutrients, temperature, etc.), have been related to MC in hundreds of studies throughout the globe and are included in the legislation of several countries (European Commission 2000; Bonada et al. 2006). In this study, we demonstrate that not only habitat and physical-chemical parameters but also pesticide residues played an important role in the composition of the macroinvertebrate community along the RMD basin (Table 6).
Most molluscs (Bivalvia, Melanoides tuberculata, Gundlachia radiata, Pyrgophorus parvulus) were associated with higher herbicide and insecticide concentrations, appearing more tolerant to such substances. The same was true for annelids (Oligochaeta and Hirudinea) and insects  (Leptoceridae, Trichoptera), and Gomphidae (Odonata), were negatively correlated with pesticide exposure (Fig. 3). These findings correspond well with the relative sensitivity to organophosphates calculated by Rico and Van den Brink (2015), except for Polycentropodidae, Caenidae, and Chironomidae, which appear very sensitive in their study (according to traits such as higher physiological sensitivity to pesticides, less migration ability, or being univoltine), but were present in highly polluted waters of the RMD. Other comparable families (listed in both studies), such as Gomphidae, Hydropsychidae, and Leptophlebiidae, were also predicted to be sensitive and correlated negatively with pesticide exposure in the present study. Information regarding higher risk traits is not available for several taxa, which could not be compared; however, molluscs, even though they were not always the same families, are predicted to be more tolerant in the study of Rico and Van den Brink (2015), which is coherent with our investigation, as well as with past research , where the gastropod Bythinella sp. appeared to be tolerant to pesticide pollution at a nearby location.
SPEAR (SPEcies At Risk) approaches (Liess and Von der Ohe 2005;Liess et al. 2008) also identified taxa that are at a higher risk of being affected by pesticide pollution. Molluscs were consistently associated with higher tolerances to pesticides, as well as annelids, the amphipod Corophiidae, and Dipterans Ceratopogonidae and Chironomidae. Leptoceridae is a SPEAR family, which was coherent with our results as well. However, main inconsistencies with respect to the present study were observed for Elmidae, Hydropsychidae, and Gomphidae, which are considered species not at risk by Liess and Von der Ohe (2005), and also for Caenidae, Polycentropodidae, and Hydroptilidae, which were  Only variables with a p value lower than 0.1 are shown. The significance levels were assessed by the Monte Carlo permutation test under the RDA option using sampling date as covariables considered species at risk, on opposition to our findings. More information with respect to traits such as life cycle (voltinism) and sensitivity to pesticides need to be gathered for tropical macroinvertebrates in order to refine predictions of their risk. Figure 4 shows the results of the variance partitioning between the three groups of variables: habitat-related parameters, physical-chemical parameters, and pesticides. Here, the highest percentage of explained variability in taxa composition results from the combination of habitat parameters and pesticides (33 %) and the combination of physical-chemical parameters and pesticides (26 %). This relates with the findings of Liess and Beketov (2011), who stated that vulnerable species with added stress were affected more strongly than non-stressed species, implying that the environmental context also determines species sensitivity to toxicants.
Arias-Andrés et al. (2016) indicated a list of high-risk pesticides for the study area. Four of those pesticides (diuron, ametryn, ethoprophos, and diazinon) are consistently highlighted as risky irrespective of the model (PERPEST, SSD, msPAF), and these are also significant explanatory variables for MC distribution in this study. It is not surprising to see very toxic insecticides (even in low concentrations) associated with diversity loss or community composition changes (Liess and Von der Ohe 2005;Bereswill et al. 2013;Ippolito et al. 2015;Rico and Van den Brink 2015); however, in this study, 50 % of the individual pesticides that correlated significantly with the MC taxa composition (Table 6) were fungicides. This should be noted, and more investigation is necessary in this area in order to calculate risk, since most fungicides act as general biocides and some can be very toxic to invertebrates (Maltby et al. 2009;Daam et al. 2009). Furthermore, initiation of downstream drift (Beketov and Liess 2008) and indirect effects such as reduced competition for food (due to death or migration of sensitive taxa) or reduced shredding activity have been associated with fungicide pollution (Rasmussen et al. 2012).
On the other hand, Rasmussen et al. (2012) discovered that streams with habitat degradation could respond with a shift in community structure that resembles that of pesticide effects, measured as the percentage of species at risk (Liess and Von Der Ohe 2005), even when there is no pesticide pollution.
Hence, the apparent response of MC to pesticide pollution in a specific site would be enhanced if there is added habitat degradation. This last statement could be the case for RMD, since degradation of habitat and some signs of organic pollution are also present within impacted sites. As usual in field studies, effects are presented as a response to multiple interacting factors, which was validated with the RDA.
Additionally, other very toxic substances, such as neonicotinoids, which are not measured in this study, could account for some of the observed variability in the MC.

Conclusions
The presence of pesticides and evidence of runoff from agricultural activities were constant at the RMD sites (Table 3, Fig. 2) during the study period.
A combination of factors related to human activities, such as habitat degradation, loss of riparian forest, altered physical and chemical parameters, and pesticide residues, is interacting in the processes that change the aquatic biological communities in the RMD.
Within the three levels of biological organization studied in this investigation, only the MC structure could be clearly related to the abovementioned stress factors (Fig. 3, Table 6). The results clearly indicate that the aquatic biota is under strong chemical pressure in the RMD watershed and that the contribution of pesticides in structuring aquatic invertebrate communities is relevant in this biodiversity-rich ecosystem.
Caño Azul (site 2) is probably the site of highest concern and one of the main contributors to the contamination of downstream ecosystems of the RMD. Based on the occurrence of high concentrations of pesticides and the toxicity of these substances in the surface waters, efforts should be assumed to assess effects in lower trophic levels since risk estimations for the RMD in companion papers submitted to ESPR (this special issue) indicate probable effects on primary producers, due to the high presence of herbicides.
Furthermore, the likelihood that the concentrations found are associated to acute and chronic toxicity events, with local species occurring in the area, is realistic (Fig. 2). Measures to mitigate pesticide runoff are crucial in order to avoid pollution of the waterways. The impact of agriculture on the environment needs to be brought into perspective in relation to the ecologic, economic, and social value of the natural habitat.
Acknowledgments The authors thank Julio Knight and his family for their valuable help during the sampling campaigns. We are grateful to the EARTH University for allowing us to take the water and biota samples inside of their campus. Also, we thank Geannina Moraga for the elaboration of Fig. 1. Jennifer Crowe and David Lean kindly revised the correct use of language. This study was funded by the Universidad Nacional de Costa Rica and the Swedish Research Council Formas, grant 205-473-3035-21.