Experimental design and community composition
We conducted a randomized-block experiment at the Russell E. Larsen Agricultural Research Center (Pennsylvania Furnace, PA, USA) with replicated mesocosm ponds. Mesocosms were 1100-L cattle tanks covered with 60% shade cloth. The spatial block was distance from a tree line in our mesocosm field. Three weeks before pesticide application, these mesocosms were filled with 800 L water, 300 g mixed hardwood leaves, and inoculations of zooplankton, periphyton, and phytoplankton homogenized from four local ponds. Just before pesticide application on the same day, each tank received two snail, three larval anuran, one larval dragonfly, one water bug, one water beetle, one larval salamander, and one backswimmer species (11 Helisoma (Planorbella) trivolvis, 10 Physa gyrina; 20 Hyla versicolor, 20 Lithobates palustris, 20 Lithobates clamitans; 2 Anax junius; 2 Belostoma flumineum; 5 Hydrochara sp.; 3 Ambystoma maculatum; 6 Nototeca undulata) (Fig. 1b). These community members naturally coexist and were applied at naturally occurring densities40. Initial conditions of some mesocosms varied in simulated pesticide treatments (see below).
We randomly assigned 18 treatments (12 pesticides, 4 simulated pesticides, 2 controls) with four replicate mesocosms of each treatment, which resulted in 72 total mesocosms (Fig. 1a). The 12 pesticide treatments were nested; we included two pesticide types (insecticide, herbicide), two classes within each pesticide type (organophosphate insecticide, carbamate insecticide, chloroacetanilide herbicide, triazine herbicide), and three different pesticides in each of four classes (Fig. 1a). To represent runoff of pesticides into freshwater systems following a rainfall event, we applied single doses of technical grade pesticides at environmentally relevant concentrations at the beginning of the experiment. To ensure our exposures represented environmental relevance, we used estimated environmental concentrations of pesticides, calculated by U.S. Environmental Protection Agency’s GENEEC v2 software, Supplementary Table 2). Our design also included water and solvent (0.0001% acetone) controls (Fig. 1a). Pesticides were obtained from ChemService (West Chester, PA, USA). Nominal concentrations of pesticides (μg/L) were: 64 chlorpyrifos, 101 malathion, 171 terbufos, 91 aldicarb, 219 carbaryl, 209 carbofuran, 123 acetochlor, 127 alachlor, 105 metolachlor, 102 atrazine, 202 simazine, and 106 propazine. We collected composite water samples 1 h after application to mesocosms and shipped samples on ice to Mississippi State Chemical Laboratory to verify these nominal concentrations. Measured concentrations of pesticides (μg/L) were: 60 chlorpyrifos, 105 malathion, 174 terbufos, 84 aldicarb, 203 carbaryl, 227 carbofuran, 139 acetochlor, 113 alachlor, 114 metolachlor, 117 atrazine, 180 simazine, and 129 propazine.
The four simulated pesticide treatments were top-down or bottom-up food web manipulations intended to mimic effects of actual herbicides and insecticides on community members. These manipulations occurred once and were concurrent with the timing of pesticide applications. Top-down and bottom-up simulated insecticide treatments were designed to reduce densities of zooplankton, simulating effects of insecticides on zooplankton survival. For top-down simulated insecticides, we doubled the densities of zooplankton predators by including six total A. maculatum larval salamanders and 12 N. undulata backswimmers per mesocosm. For bottom-up simulated insecticides (i.e., direct manipulation of a lower arthropod trophic level), we removed zooplankton with a net. Top-down and bottom-up simulated herbicides were designed to reduce algae, simulating effects of herbicides on survival and growth of algae. For top-down simulated herbicides, we doubled the densities of large herbivores to increase grazing pressure by including 22 H. trivolvis snails, 20 P. gyrina snails, 40 H. versicolor larval anurans, 40 L. palustris larval anurans, and 40 L. clamitans larval anurans per mesocosm. For bottom-up simulated herbicides, we covered mesocosms in three sheets of 60% shade cloth in an attempt to block light and reduce photosynthesis. The experiment ran for four weeks, from June to July.
Measurements of experimental responses
During the experiment, we sampled periphyton using clay tiles (100 cm2) oriented perpendicularly along the bottom of the mesocosm. Each mesocosm had two periphyton measurements: ‘inaccessible periphyton’ taken from caged clay tiles that excluded herbivores and ‘accessible periphyton’ taken from clay tiles that were uncaged, allowing herbivore access. We sampled phytoplankton from water samples taken 10 cm below the water surface. Periphyton was scrubbed from tiles and phytoplankton from water samples (10 mL) were filtered onto glass fiber filters (under low vacuum pressure, <10 psi; Whatman EPM 2000, 0.3 μm, 47 mm) to estimate associated chlorophyll concentrations. The chlorophyll concentration of each filter was determined using an organic extraction procedure with a 50:50 mixture of 90% acetone to DMSO. We measured chlorophyll-a concentrations using a standard fluorometric technique. We scored water clarity, a metric of light availability, on a scale from one (clear) to five (opaque) blinded to treatment. We measured pH and dissolved oxygen (DO) at dusk and dawn on subsequent days using hand-held meters (YSI, Yellow Springs, OH, USA). We measured decomposition by taking the dry mass of hardwood leaf packets in each mesocosm at the beginning and the end of experiment. In addition, we sampled snail egg masses and hatchlings using two rectangular pieces of Plexiglass (465 cm2) in each mesocosm, one hung on the side and one on the bottom of the mesocosm. Zooplankton were collected from the entire water column by placing a PVC pipe (10 cm diameter, 60 cm height) upright in the center of each tank, capping the bottom, and pouring the water through a 20 μm Nitex mesh. We collected two samples of zooplankton from each mesocosm, and we combined and preserved the samples in 70% ethanol. Zooplankton were counted and identified in 5 mL subsamples for each mesocosm using a zooplankton counting wheel (Wildlife Supply Company, Yulee, FL, USA) and a dissecting microscope. At the end of the experiment, mesocosms were drained, and the remaining animals were counted, euthanized, and preserved. Two previous manuscripts, which use the same design as the current manuscript, also describe this experimental design and methods8,41. The research was reviewed and approved by the Penn State University Institutional Animal Care and Use Committee.
Statistical analyses
To test for the consistency of effects of type, class, and individual pesticide on aquatic ecosystem processes and communities and to attribute the variation explained to each pesticide level of organization while accounting for the nested structure of our experimental design (Fig. 1a), we completed permutational analyses of variance (PERMANOVA). For nested PERMANOVA models, the predictors were the following random categorical terms: type (insecticide, herbicide), class (carbamate, organophosphate, chloroacetanilide, triazine) nested within type, and pesticide (12 in total) nested within class within type. These models did not include controls or simulated pesticides because these treatments were not hierarchically nested (Fig. 1a). We evaluated 9999 permutations using residuals under a reduced model. Following nested PERMANOVAs, we used pair-wise multiple comparisons tests using PERMANOVAs to evaluate differences among controls, organophosphates, carbamates, top-down simulated insecticides, bottom-up simulated insecticides, chloroacetanilides, triazines, top-down simulated herbicides, and bottom-up simulated herbicides. In these pair-wise comparisons, we evaluated 9999 unrestricted permutations of raw data. All PERMANOVAs also included spatial block as a random predictor to account for variation in sunlight associated with distance from a tree line. Preliminary analyses showed that exclusion of the block did not change the results. In all PERMANOVAs, test statistics associated with Type III partial sums of squares were evaluated.
We conducted four nested PERMANOVAs. Our first nested PERMANOVA focused on ecosystem processes and included the following responses: pH taken at dawn, respiration (the difference between dissolved oxygen at dusk and dissolved oxygen at dawn of the subsequent day), decomposition (percent mass remaining of hardwood leaf packets), turbidity (water clarity scores from 1 to 5), and densities of phytoplankton and accessible periphyton (measured via chlorophyll-a). Since this analysis was focused on ecosystem-level responses, for periphyton we include accessible periphyton and not inaccessible periphyton because accessible periphyton better encompasses the totality of biologic (e.g. herbivore predation, shading from phytoplankton) factors that could influence periphyton biomass at the ecosystem level. Preliminary analyses included both accessible and inaccessible periphyton, and the results did not differ. The resemblance matrix for these responses was constructed using a Euclidean distance matrix of log-transformed and normalized values.
Our second and third nested PERMANOVAs focused on community structure. We separated community members into two statistical models based on the forms of response variables; those whose response variables were densities based on counts (zooplankton) and those community members whose response variables were survival, mass, reproductive rates, or density abstracted from chlorophyll measurements (insect predators, snail and tadpole herbivores, and algae; termed the tri-trophic community). The multivariate response for the zooplankton community included densities of Daphnia, Diaphanasoma, Chydorus, Bosmina, Diaptomus, and Cyclops. Zooplankton community analyses were based on square-root transformed densities using Bray-Curtis similarities. The multivariate response for the tri-trophic community model included: survival (0 to 1) of all amphibian, snail, and insect community members; average masses of surviving individuals for each amphibian species and H. trivolvis snails; average number of hatchlings and eggs per surviving H. trivolvis snail; and densities of phytoplankton and accessible and inaccessible periphyton. In this analysis, we include both accessible and inaccessible periphyton to account for any differential effects of excluding herbivores on the abundance of periphyton. Mass and reproductive rates were standardized to the number of surviving individuals to account for the different densities added to each tank at the beginning of the study (i.e. extra herbivores in top-down simulated herbicide treatment and extra predators in bottom-up simulated insecticide treatment). Finally, our fourth nested PERMANOVA evaluated a simplified tri-trophic community. We simplified the tri-trophic community responses into three functional roles within the community: algae, herbivores, and predators. Tri-trophic community responses of individua taxa were transformed and normalized as described previously, and then they were averaged according to functional group. We averaged densities of periphyton and phytoplankton into a single “algae” response, all amphibian and snail responses into a single “herbivore” response, and all insect and salamander responses into a single “predator” response. The simplified tri-trophic community model was based on Euclidean distances.
To visualize consistency of effects within type, class, and individual pesticides on multivariate ecosystem and community responses and to compare pesticide effects to simulated pesticides and controls, we used distance-based redundancy analyses (dbRDA) and two-way cluster diagrams. The dbRDAs were based on appropriate resemblance matrices for ecosystem and community responses as described above. The underlying categorical predictors in all models included: the spatial block, organophosphate, carbamate, chloroacetanilide, triazine, top-down simulated insecticide, bottom-up simulated insecticide, top-down simulated herbicide, bottom-up simulated herbicide, and control. In the dbRDA plots, we show the centroid values for the 18 experimental replicates.
As an alternative to the dbRDAs presented in the main text, we also visualized the consistency of effects within type, class, and individual pesticide on ecosystem, tri-trophic community, and zooplankton responses and compared pesticide effects to simulated pesticides and controls, using principal coordinates analyses (PCoA) (Supplementary Figs. 1, 2, 4). PCoAs were based on appropriate resemblance matrices as described previously. PCoAs were conducted in PERMANOVA+ for PRIMER and resulting data were exported. Point and vector plots were made using the exported data and the ‘ggplot2’ package in R. Ellipses on point plots represent 95% confidence intervals of groups based on standard errors and were made using the ordiellipse function in the ‘vegan’ package.
For the two-way cluster diagrams, clusters of pesticide treatments were based on centroid distances of the appropriate resemblance matrices. Clusters of multivariate responses were based on Euclidean distance resemblance matrices of averaged treatment responses. Before averaging, ecosystem and community responses were transformed and normalized as described previously. In clustering of treatments and responses, the cluster mode was the group average. In the PERMANOVAs for tri-trophic community responses, the effect of block was significant (Supplementary Table 1). Thus, we accounted for the effect of block by taking the residuals of simple linear regressions with individual tri-trophic community responses as the independent variable and block as the predictor in the generation of the shaded values of the two-way cluster diagrams. Then, we averaged these block-adjusted treatment responses with the ‘shade plot’ function in PRIMER. For the ecosystem responses and zooplankton community, the effect of block was not significant in the PERMANOVA model (Supplementary Table 1), so shaded values of the two-way cluster diagrams were simply the averaged treatment responses. All PERMANOVA models, pair-wise comparisons, dbRDAs, and two-way cluster diagrams were executed using PERMANOVA+ for PRIMER version 7 (PRIMER-E Ltd, Plymouth, UK). For ease of visualization of dbRDA and PCoA plots, data from PERMANOVA+ for PRIMER were exported, and plots were made using ‘ggplot2’ package in R.
To quantitatively test whether changes in composition, abundance, and/or richness of different functional groups mediate the effects of herbicides and insecticides on ecosystem function, we performed structural equation modeling and model comparison using the ‘piecewiseSEM’ package. We conducted separate analyses for insecticides and herbicides. Global models, which replicated the structures of periphyton and phytoplankton food webs (Fig. 1b), evaluated how herbicides or insecticides influence ecosystem functions through changes in composition, abundance, and richness of: (1) snail and tadpoles predators, (2) snails and tadpoles, (3) zooplankton predators, and (4) zooplankton. Ecosystem functions included primary production of accessible periphyton, primary production of phytoplankton, and whole system respiration. We did not include decomposition because ordination analyses showed decomposition was not dramatically influenced by pesticide exposures. Since this analysis was focused on ecosystem-level responses, we included accessible periphyton only (i.e. excluded inaccessible periphyton) because accessible periphyton better encompasses the totality of biologic (e.g. herbivore predation, shading from phytoplankton) factors that could influence periphyton biomass at the ecosystem level. For each functional group, composition was the locations of species or genera along the first axis of a principal coordinates analysis (PCoA) using a Bray-Curtis similarity matrix. Abundance was the total number of individuals of a functional group in a mesocosm, and richness was the number of species or genera within a functional group. The relationships among these variables in the global models are presented in Supplementary Fig. 5. Notably for herbicides, the global model considered the bottom-up direct effects of herbicides on periphyton and phytoplankton and the non-target direct effects of herbicides on zooplankton composition, abundance, and richness. Mesocosms treated with chloroacetanilide herbicides were excluded from the analysis because they did not influence algal densities (Fig. 3c). For insecticides, the global model considered the top-down direct effects of insecticides on composition, abundance, and richness of tadpole and snail predators, zooplankton predators, and zooplankton. These global models did not consider non-target direct effects of insecticides and herbicides on tadpoles and snails because community tri-trophic community analyses (Fig. 3) showed little evidence of non-target effects on these organisms. Phytoplankton, periphyton, and abundance of functional groups were all log-transformed in models. To facilitate comparisons among responses and clarify relationships among predictors, we reduced the global models by removing non-significant paths. The final models are presented in Fig. 4. All data supporting the results of the present manuscript are publicly available42. Full information regarding creation and attribution of silhouettes used throughout the main text and supplementary figures is included in Supplementary Table 3.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
