Effects of multiple stressors associated with agriculture on stream macroinvertebrate communities in a tropical catchment

Tropical forests are declining at unprecedented rates in favour of agriculture, and streams can be severely impacted due to effects of multiple stressors that have rarely been considered together in tropical studies. We studied the effects of multiple stressors associated with agricultural practices (pesticide toxicity, nutrient enrichment and habitat alteration–quantified as TUmax, soluble reactive phosphorus concentration and sedimentation, respectively) on macroinvertebrate communities in a tropical catchment in Panama (13 stream sites sampled in 20 occasions from 2015 to 2017, with 260 samples in total). We examined how macroinvertebrate abundance, taxonomic richness, community composition and biotic indices (SPEAR and BMWP/PAN, which were specifically designed to detect pesticide toxicity and nutrient enrichment, respectively) varied depending on the studied stressors, considering their single and combined effects. Our analyses revealed significant effects of the studied stressors on macroinvertebrate communities, with two particular results that merit further attention: (1) the fact that pesticide toxicity affected BMWP/PAN, but not SPEAR, possibly because the former had been adapted for local fauna; and (2) that most stressors showed antagonistic interactions (i.e., lower combined effects than expected from their individual effects). These results highlight the need for toxicity bioassays with tropical species that allow adaptations of biotic indices, and of observational and manipulative studies exploring the combined effects of multiple stressors on tropical macroinvertebrate communities and ecosystems, in order to predict and manage future anthropogenic impacts on tropical streams.

Introduction on average and up to 7,000 mm on the highlands, with 87.7% occurring in the rainy season (May-December) [24].
The study catchment is intensely used for agriculture, being one of the most productive areas in Panama [26]. The strong erosion, as a result of native vegetation removal, steep slopes and high precipitation, causes the progressive deterioration of the catchment, and stream water quality is affected by the entrance of fine sediment, pesticides and nutrients, the latter coming both from fertilizers and from the inefficient treatment of waste water in the area [27]. We conducted the study at 13 sites (Fig 1; S1 Table)  mm), and water chemistry; (3) collected water samples for further physico-chemical analyses and determination of pesticides; and (4) sampled macroinvertebrates.
We characterized the habitat using the rapid habitat assessment protocol developed by Barbour et al., [28] for the United States Environmental Protection Agency (EPA) for high-gradient streams. This consisted of qualifying 10 variables (epifaunal substrate/available cover, embeddedness, velocity/depth regime, sedimentation, channel flow status, channel alteration, frequency of riffles, bank stability, bank vegetative protection, and riparian vegetative zone width) using a numerical scale from 0 to 20 (maximum). Each variable was assessed independently, and their sum was assigned to one of four categories of habitat quality (i.e., optimal, suboptimal, marginal or poor).
Substrate composition was characterized visually as the proportion of different size classes of mineral substrate (boulder, cobble, gravel, coarse and fine sand, and clay) and CPOM and FPOM were quantified visually as the proportion of streambed covered by each type of organic matter [28]. We measured pH, temperature (˚C), conductivity (μS cm -1 ), turbidity (mg L -1 ) and dissolved oxygen (% saturation) in situ using a multiparametric probe (YSI 556), and collected two sets of 2-Lwater samples from the mid column in middle of the stream, which were transported to the laboratory on ice. The first set of water samples was analysed at the Environmental Quality Laboratory of the Ministry of Environment (Panama) for concentrations (mg L -1 ) of total solids using a gravimetric method (SM 2540 B), and nitrate and soluble reactive phosphate (SRP) using spectrophotometric methods (SM 4500-NO 3 B and SM 4500-P B5 and E) [29].

Determination of pesticides
The second set of water samples was analysed for pesticides at the Plant Health Laboratory from the Agricultural Development Ministry (MIDA, Panama). A 2-L water sample was collected at each site from the middle of the stream and the mid column. Samples were immediately refrigerated and transported to the laboratory, and kept at 4˚C until analysis was performed within 24 h of receipt. Pesticides were determined using two methods: liquid-liquid microextraction [30] and direct injection [31]. The first method was used for organophosphates, organochlorines and pyrethroids; pesticides were extracted with ethyl acetate and residuals were quantified by gas chromatography and mass spectrophotometry (GC-MSMS; limit of quantification: 0.11 μg L -1 ). The second method was used for triazines, carbamates and other polar pesticides; samples were injected and analysed with high performance liquid chromatography and mass spectrophotometry (LC-MSMS; limit of quantification: 0.10 μg L -1 ) and electrospray ionization with dynamic acquisition (MRM mode), which avoids solid phase extraction. The percentage of recovery ranged between 70 and 110% (CV = 11%). Linearity was measured by the R 2 coefficient for the individual pesticide calibration curves, always resulting in R 2 � 0.99. Each set of samples was analyzed in duplicate, simultaneously with a laboratory blank. To avoid matrix effects we used a matrix-matched calibration curve.

Macroinvertebrate sampling and processing
Macroinvertebrates were kick sampled using a 30-cm wide D-net with a 0.5-mm mesh. At each site we took three 2-m long samples, which were subsequently pooled, with a total area of 1.8 m 2 sampled per site. Samples were taken on a variety of habitats including mineral substrate, leaf litter patches and bank vegetation, in proportions similar to their presence in the stream. The net contents were transferred to a 0.5-mm mesh sieve and then to a white tray, where macroinvertebrates were preliminary sorted, and stones, leaves and wood discarded. The rest of the sample was introduced in labelled vials with 96% ethanol and transferred to the Freshwater Macroinvertebrate Laboratory at the COZEM-ICGES (Panama). Macroinvertebrates were sorted and identified to family level-which is the usual procedure to calculate the SPEAR and BMWP indices [18,32]-using identification keys for tropical taxa [33][34][35][36].

Calculation of pesticide toxicity
In order to have a standard value of toxicity associated with pesticide concentrations measured at each site we used the Toxic Unit (TU) approach [37]. The TUs were given as maximum TU (TU max ), a simple approach widely used in the literature [15;38,39]. To calculate TU max we considered all pesticides found across samples at each site, excluding those below the quantification limit. Given that toxicity data for tropical stream macroinvertebrates are unavailable, we calculated TU max based on data available for Daphnia magna [15] based on the following equation: where TU (D. magna) is the TU max of n pesticides detected in the study site, C i is the concentration of pesticide i (μg L -1 ), and LC50 i is the 48 h acute median lethal concentration (μg L -1 ) reported for pesticide i in D. magna.

Calculation of biotic indices
To calculate the SPEAR index, taxa were classified into species at risk (SPEAR) or species not at risk (SPEnotAR) according to several ecological and physiological traits [15], which were obtained from an online database (http://www.systemecology.eu/spear/spear-calculator/). The SPEAR value for each site was calculated as follows: where n is the number of taxa, x i is the abundance of taxon i, and y is 1 if taxon i is classified as SPEAR, otherwise 0. The BMWP index is one of the most often used indices based on macroinvertebrates to assess nutrient enrichment in streams [40]. It was originally developed for the United Kingdom [18] and has been adapted to the local fauna of many countries, including Panama (BMWP/PAN; [27]). The BMWP/PAN was adapted based on tolerance to nutrient enrichment of local macroinvertebrate families, following the methods of Ruiz-Picos et al., [40]. The BMWP score at a given site is the sum of the individual scores of the families present at that site, which range from 1 (most tolerant families) to 9 (most sensitive families).

Statistical analyses
All analyses were performed in R software, version 3.6.0 [41]. We first explored bivariate scatterplots and Pearson correlations to select the most relevant and uncorrelated environmental variables (r � 0.70) to be used in further analyses (S1 Fig; [ 42]); these variables were TU max (hereafter pesticide toxicity), SRP concentration (hereafter nutrient enrichment), the sediment deposition index (hereafter sedimentation index; inversely related to sedimentation), and water temperature (hereafter warming); other variables were discarded. Scatterplots and correlations were performed with the "chart.Correlation" function in PerformanceAnalytics package [43].
Secondly, we examined the individual and interactive effects of pesticide toxicity, nutrient enrichment and sedimentation index on macroinvertebrate abundance, taxonomic richness and biotic indices (SPEAR and BMWP/PAN), using linear mixed-effects models accounting for temporal autocorrelation. Warming influence was not considered in these models to avoid the complexity of a four-way interaction model, and because sedimentation was a better representation of habitat alteration (S2 Table). Models were first defined in terms of random structure, and a model selection procedure was used to identify interactions between predictors [42]. The optimal model random structure (i.e., the need for a variance structure, temporal correlation structure and/or random term) was defined by comparing models containing different terms using the Akaike Information Criterion corrected for sample size (AICc) (S3 Table). Final models were fit using the "lme" function (linear mixed effects), with site as a random term (except for the richness model, which lacked this component), temporal autocorrelation (ARMA correlation structure), and a variance structure (VarIdent in relation to site to control for different variances within sites).
Interactive effects were explored through five models, all containing the three predictors, but varying in the number of interactions. The null model (model 1) assumes no interactions between predictors (i.e., additive effects only); three models (models 2, 3 and 4) included pairwise interactions between nutrient enrichment and sedimentation index, pesticide toxicity and sedimentation index, or pesticide toxicity and nutrient enrichment; and one model (model 5) included all interactions, including the three-way interaction. The five models were compared using an AICc-based model selection approach, with the most plausible models being selected based on delta AICc (Δi; i.e., difference in AICc value relative to the best model) and Akaike weights (wi; i.e., the probability that a model is the best among the whole set of models). Residuals from each model were inspected to ensure there were no visual patterns and that linear model assumptions (i.e., independence and homogeneity assumptions) were not violated. Estimates and 95% confidence intervals for single predictors and their interactions were obtained using a model averaging approach, which averages the estimates of the retained models containing the parameter. Models were constructed, selected and averaged using nlme ("gls", "lme", "VarIdent" and "corARMA" functions; [44]) and MuMIn ("model. sel" and "model.avg" functions; [45]) packages.
Thirdly, we evaluated the effect of pesticide toxicity, nutrient enrichment and habitat alteration (sedimentation index and warming) on macroinvertebrate community composition using redundancy analysis (RDA; [46]), where the species dataset was predicted by the environmental dataset. Both datasets contained multiple samples taken over time and were averaged to produce a single value per site. Lastly, to quantify the amount of variability in macroinvertebrate community composition that can be attributed to each of the above environmental factors, as well as to their shared contribution (i.e., interactions between predictors), we used partial redundancy analysis [47]. The amount of variability explained by each factor and their shared contribution was based on adjusted R 2 (R 2 adj ), and their statistical significance tested through permutation tests (999 randomizations). Macroinvertebrate data were Hellinger-transformed prior to both procedures to provide an unbiased estimate of variance partitioning based on pRDA. Variance partitioning and permutation tests were performed using the "varpart" and "cca.anova" functions, respectively, both from the vegan package [48]. Results were presented using a Venn diagram, which was drawn on Inkscape, an open-source vector graphics editor.

Interactive effects of pesticide toxicity, nutrient enrichment and habitat alteration on macroinvertebrate communities
The model selection procedure revealed that, in most cases, there were two best models explaining the observed patterns (~60% probability based on Akaike weights); the exception was abundance, which was explained by a single model with pairwise interactions. The SPEAR and BMWP indices were best explained by one additive model (i.e., without interactions) and one model containing pairwise interactions; the two most plausible models explaining richness contained pairwise interactions (Table 1; S9 Table). Overall, individual effects of pesticide toxicity and nutrient enrichment were negative, while the sedimentation index had a positive effect (indicating a negative effect of sedimentation, which was inverse to the index). The sedimentation index was the only factor individually affecting all the response variables; the individual effect of nutrient enrichment was important for both biotic indices, but not for abundance or richness; and pesticide toxicity individually affected all variables except SPEAR (Fig 2). The interaction between pesticide toxicity and sedimentation index was significant for abundance, richness and BMWP, always having a negative antagonistic effect (i.e., lower than predicted by the sum of individual effects); the interaction between nutrient enrichment and sedimentation index was important for richness and SPEAR, with a positive additive and a negative antagonistic effect, respectively (Fig 2).
All stressors explained 62% of variance in macroinvertebrate community composition. Nutrient enrichment and sedimentation were mostly related to RDA1 (both with positive relationships; the sedimentation index being inversely related to sedimentation), while pesticide toxicity and warming were related to RDA2 (both with negative relationships) (Fig 3). Thus, some taxa were related to sites with lower levels of pesticide toxicity, nutrient enrichment and habitat alteration (i.e., sites S-02, S-04 and S-10; Hyalellidae, Leptophyphidae, Leptophlebiidae, Planariidae, Planorbidae, Ptilodactylidae, Odontoceridae and Tabanidae) and others were associated to more impacted sites, that is, affected by nutrient enrichment and sedimentation (i.e., S-08 and S-12; Baetidae and Hydroptilidae) or higher levels of pesticide toxicity and warming (i.e., S-06 and S-07; Chironomidae, Lumbriculidae and Psychodidae).
The pRDA showed that a large proportion of variance in macroinvertebrate communities was driven by nutrient enrichment (R 2 adj = 0.51) and habitat alteration (R 2 adj = 0.37), while pesticide toxicity contributed to a lower proportion of variance (R 2 adj = 0.13). The proportion of variance attributable to the combination of pesticide toxicity and nutrient enrichment (R 2 adj = 0.50), nutrient enrichment and habitat alteration (R 2 adj = 0.46) or the whole set of environmental predictors (R 2 adj = 0.62), was lower than expected based on the sum of individual stressor effects (i.e., the additive expectation), indicating antagonistic effects. On the other hand, the combination of pesticide toxicity and habitat alteration (R 2 adj = 0.33) was slightly higher than expected, suggesting a synergism between these two stressors ( Table 2; Fig 4).

Discussion
Assessing the effects of agricultural practices on tropical stream communities is an urgent challenge, given the fast conversion of tropical forests to agricultural land due to the rising demands of human populations [4,5,49]. Studies, however, are scarce and have only partially addressed this question, as they have focused on single stressors such as pesticide toxicity Table 1

. Summary of model selection testing for interactions between multiple stressors on macroinvertebrate abundance and richness and the SPEAR and BMWP indices, based on the Akaike information criterion corrected for sample size (AICc).
Models are ordered from the best to the poorest fit according to Akaike weights (wi). K, number of estimated parameters for each model; Δi (delta AICc), difference in AICc value relative to the best model; wi, probability that a model is the best among the whole set of models. For each response variable, five models were constructed, which are ordered from the simplest model without interactions (model 1: null model, with no interactions) to the most complex one (model 5, containing the 3-way interaction). Models differ in the number of parameters according to the most parsimonious combination of structure and terms described in S9 Table. PT, pesticide toxicity (Tu max ); SE, sedimentation index; NE, nutrient enrichment (SRP).  [6,50,51], nutrient enrichment [52] or habitat alteration, mainly deforestation [53][54][55] and sedimentation [20,56]. Our study is, to our knowledge, the first to assess the joint effect of multiple stressors associated with agriculture on tropical stream macroinvertebrate communities. We demonstrated negative effects of the studied stressors (pesticide toxicity, nutrient enrichment, sedimentation and warming) on macroinvertebrate community descriptors and/ or biotic indices. Sedimentation was the only factor with negative effects on all the variables; this factor has been shown to have large effects on tropical macroinvertebrates, which move downstream in response to increased sedimentation [20]. Abundance and richness were not affected by nutrient enrichment, in agreement with other tropical studies and possibly because other factors (e.g. light) limited primary productivity [12]. In contrast, abundance and richness were negatively affected by pesticide toxicity, an effect that has not been found elsewhere in the tropics [6,50]. The different stressors caused shifts in community composition, with some taxa being more tolerant to pesticide toxicity or warming (i.e., some dipterans and oligochaetes) and others to nutrient enrichment or sedimentation (i.e., some mayflies and caddisflies).

Model
Interestingly, pesticide toxicity affected the BMWP/PAN but not the SPEAR index, which had been specifically designed to assess pesticide effects on macroinvertebrates [15]. This may   Table 2. Results of partial redundancy analysis (pRDA). Exploring the amount of variance in macroinvertebrate community composition explained by pesticide toxicity (TU max ), nutrient enrichment (SRP) and habitat alteration (temperature and sedimentation index). We shown the degrees of freedom (df model , df residual ), adjusted R 2 (R 2 adj ), associated p-values (p; after permutation tests using 999 randomizations), additive expectation (sum of R 2 adj of individual stressors), and interaction type (A; antagonistic when R 2 adj of interaction is lower than the sum of individual stressors; S, synergistic when R 2 adj of interaction surpasses the sum of individual stressors be due to the fact that the SPEAR index is based on physiological traits associated with pesticide sensitivity in temperate species, which highlights the need for conducting biological toxicity tests with tropical macroinvertebrates, as these are likely to show different environmental sensitivities even at the taxonomic resolution of family [57]. This is supported by the fact that only studies in temperate areas have shown reduced levels of SPEAR with increased pesticide toxicity [15,58,59]. In our study, BMWP/PAN was affected by all the studied stressors, including pesticide toxicity. While the BMWP index was designed to assess effects of nutrient enrichment on macroinvertebrates [60], we used an index that had been adapted for local fauna (in contrast to SPEAR) and statistically calibrated [27,40], which may explain its significant response to all stressors. Temperate studies have also found an effect of pesticide toxicity on BMWP (but see [16,61]), while this has not been the case for other tropical studies using adapted versions of the index, such as the BMWP/COL [51].
Importantly, our analyses revealed interactive effects of different stressors that, in most cases, were antagonistic. Specifically, effects of pesticide toxicity or nutrient enrichment in combination with sedimentation on community descriptors and biotic indices were lower than expected from single effects, and the combined effects of most stressors on community composition were also antagonistic. These results suggest that assessing effects of stressors associated to agriculture individually can overestimate their overall effect, and highlights the importance of using a multi-stressor approach in real-context studies, because of the complex and often unpredictable interactions between stressors [10]. Our results are in accordance with a recent literature review, which found that additive and antagonistic interactions of multiple stressors were more prevalent than synergistic interactions [62]. Further studies should explore interactions of co-occurring stressors in the field, but also under controlled conditions where stressors can be easily manipulated (e.g., fully factorial designs in microcosms or mesocosms).
In summary, we provided novel evidence about negative effects of agricultural practices on tropical stream macroinvertebrate communities, which were affected by multiple stressors acting in combination. Our results highlight the need for (1) further tropical studies using a multi-stressor approach, including observational and manipulative studies assessing how macroinvertebrate communities and ecosystems respond to different combinations of stressors and; and (2) toxicity bioassays with tropical species that allow the adaptation of biotic indices to local fauna. Moreover, functional metrics such as leaf litter breakdown or other ecosystem processes can be useful tools for detecting ecosystem responses to nutrient enrichment [63,64], although these metrics also respond to other stressors and environmental drivers. Thus, the combined use of structural and functional metrics (e.g., biotic indices and ecosystem-level processes) could provide a more comprehensive assessment of the ecological effects of multiple stressors [65]. This is particularly needed in tropical areas, which are understudied and subject to rapid transformation by human activities [5], and whose responses compared to their temperate counterparts are difficult to predict [66].

S1 Fig. Correlogram showing pairwise Pearson correlations between variables sampled.
Positive correlations are displayed in blue and negative correlations in red color. Color intensity and the size of the circle are proportional to the correlation coefficients. In the right side of the correlogram, the legend color shows the correlation coefficients and the corresponding colors. Only significant (p-value < 0.05) are displayed (blank spaces indicate non-significant correlations). (DOCX) S1  Table. Results of model selection to define the best random structure of models, in terms of variance structure, temporal correlation and random term.