Set against a backdrop of growing human-populations, most habitats are affected by more than one anthropogenic pressure (Halpern et al., 2008; Barnosky et al., 2011; Bennett et al., 2015; Halpern et al., 2015; Tilman et al., 2017) and are rarely subject to the actions of single stressors/disturbances (Folt et al., 1999; Crain et al., 2008). Research to understand the effects of multiple stressors on biological and ecological processes has been increasing in recent years (Underwood, 1980; Vinebrooke et al., 2004; Blake and Duffy, 2010; Darling et al., 2013; MacDougall et al., 2013; Côté et al., 2016; Galic et al., 2018; Mayer-Pinto et al., 2020; Jackson et al., 2021). At its most simple, multiple stressors can cause chemical, biological and ecological impacts in three ways depending on whether the combined effect of both stressors is smaller (antagonistic interaction), equal (additive) or greater than (synergistic interaction) the sum of the effects of the individual stressors (Folt et al., 1999; Crain et al., 2008; Lyons et al., 2015). Despite this simplicity, there is still debate about the terminology used to describe combined effects, how prolific each non-additive effect is and whether their frequency and magnitude modify their effects (Piggott et al., 2015; Côté et al., 2016; Thompson et al., 2018; Orr et al., 2020).

In part, this is because it is difficult to predict the effects of multiple stressors because those effects can vary depending on the intensity of each stressor involved, with some synergistic effects at smaller concentrations and others antagonistic at larger concentrations (Chen et al., 2004; Coors and De Meester, 2008; Townsend et al., 2008). Yet, little of that work has, however, considered more than one or two concentrations/intensities of each stressor with research mainly having been performed in laboratory settings with only 10-20% of studies being undertaken in the field with natural-occurring assemblages (Crain et al., 2008; Mayer-Pinto et al., 2020). In addition, field experiments with the ability to test multiple pairwise comparisons have been rare (Côté et al., 2016) despite calls for their implementation (Strain et al., 2014; Bennett et al., 2015; Lyons et al., 2015; Côté et al., 2016; Orr et al., 2020; Rindi et al., 2022). This is predominantly because the large complex experiments required have been very difficult to run in the field. A recent review on the sub-lethal effects of multiple contaminants on habitat-forming species, revealed ~40% of publications were poorly designed with insufficient replication and/or appropriate controls (Mayer-Pinto et al., 2020). A suitable experimental design for such research is the ‘response-surface’ design (Inouye, 2001; Inouye, 2005; Rosenblatt and Schmitz, 2016). ‘Response-surface’ designs test the relationship between the response and a number of input (predictor) variables across a range of levels (Box and Wilson, 1951; Inouye, 2005; Lin and Peterson, 2006; Varaksin et al., 2018). Response-surface experiments have tended to be done in the laboratory with individual species or simplified assemblages, however, they have been criticized for lacking replicate combinations of treatments (Benedetti-Cecchi, 2005) and environmental relevance (Mayer-Pinto et al., 2010). This makes it difficult to examine how chemical and biological effects on organisms can be translated into ecological effects on assemblages, which are necessary to understand the effects of multiple stressors on the ecosystem (Attrill and Depledge, 1997). Wildlife exposed to multiple stressors, and not examined using appropriate experiments, surveys and modelling can fail to identify and manage declines in populations (Munns, 2006; Browne et al., 2015; Côté et al., 2016; Rosenblatt and Schmitz, 2016). Therefore, if we are to improve legislation and guidelines so they protect habitats and the populations of species they contain, as advocated in the literature recently (Orr et al., 2020; Jackson et al., 2021), appropriately designed manipulative field experiments are needed to test the effects of multiple stressors and their interactions, across a range of intensities/concentrations and across different levels of organization within the ecosystem. Here we present a ‘response-surface’ field with replicated combinations of treatments to (a) test whether the responses of mussels (copper accumulation, cellular viability) and their assemblages (filtration, respiration) interact additively, synergistically, or antagonistically to variations in the concentration of two stressors. We hypothesize that the combined influence of stressors (copper and the biocide chlorpyrifos) will vary depending on their relative concentrations, with maximal combined influences not necessarily being predictable based on individual dose-response curves and (b) understand if individual and combined effects of multiple stressors at a range of concentrations are consistent among chemical, biological, and ecological responses derived from different levels of biological organization.

Predicting the combined effects of stressors based on knowledge of their individual effects is fraught with uncertainty; particularly when more recent knowledge suggests that non-additive effects may outweigh additive effects (Strain et al., 2014; Piggott et al., 2015; Côté et al., 2016). Using a ‘response-surface’ field experiment this study reveals a prevalence of non-additive effects between stressors that varied depending on their concentration, and on the responses of the mussels and assemblages measured. For the cellular level response variable, there were antagonistic interactions between increasing concentrations of both stressors that peaked at intermediate concentrations of chlorpyrifos and the highest concentration of copper such that the cellular viability of mussels at those concentrations were similar to the controls ( Figures 1A, B ). In terms of copper in the tissue of mussels, there was a synergistic interaction in which mussels exposed to large concentrations of copper and intermediate concentrations of chlorpyrifos accumulated significantly more copper than mussels exposed to the highest concentration of copper only ( Figures 2A, B ). Additionally for filtration, there was evidence of two non-additive effects with filtering by mussel assemblages exposed from intermediate to larger concentrations of chlorpyrifos across the gradient of copper concentrations being reduced to the same extent as the single stressor with the greatest effect; in this instance the intermediate to highest concentrations of chlorpyrifos ( Figures 4A, B ). This particular type of non-additive effect is suggestive of one that is predicted under the “simple comparative/dominance model” where a single stressor dominates the other stressor in determining their combined effect (Folt et al., 1999; Côté et al., 2016). Additionally, there was a synergistic interaction between small to intermediate concentrations of chlorpyrifos and the highest concentration of the copper, which reduced filtering by assemblages when compared to the controls ( Figures 4A, B ). In terms of respiration by assemblages, there was a significant interaction between copper and chlorpyrifos ( Table 1C ; Figures 3A, B ), however, given the small amount of deviance explained for the model (9%) and considering that p-values in a GAM are only approximate, caution must be taken in the interpretation.

One possible explanation for the antagonistic interaction between stressors occurring for cellular viability of mussels could be due to exogenous (direct) interactions between the stressors in the environment (i.e. in the water-column before entering the organism) (Lyons et al., 2015). These types of chemical interactions are thought to occur in aquatic environments between these contaminants, where chlorpyrifos may chelate copper-ions (Meikle and Youngson, 1978). Others, however, believe that copper’s speciation is not only indicative of its toxicity but also indicative of its bioavailability (Hall and Anderson, 1999). Although some more benign species of copper can still be assimilated into the tissues of bivalves through feeding or respiration, this is thought to be at much smaller concentrations than species of copper that are thought to be more toxic (Pan and Wang, 2009). As such, if there was an exogenous interaction between the stressors, one would have expected to observe a reduced uptake of both stressors when combined. While we were unable to quantify chlorpyrifos in the tissues of mussels, there was a synergistic interaction between intermediate chlorpyrifos-concentrations and increasing copper-concentrations, which resulted in mussels at those mixtures accumulating more copper in their tissue than mussels exposed to the highest concentration of copper alone ( Figure 2A ).

In fact, Dondero et al. (2011) reported similar antagonistic interactions between increasing concentrations of this chlorpyrifos and another metal (nickel) on cellular responses in mussels; however, the authors found no evidence of exogenous chemical reactions. Instead, the authors found that the antagonism was due to the stress-response of the species. They used molecular techniques to identify whether there were more genes associated with metals (metallothioneins). This latter finding may be caused by an endogenous antagonistic interaction between the contaminants, where an organism responds to different contaminants through a general response (i.e. inflammation; Browne et al., 2015) that may include the rapid synthesis of proteins (e.g. enzymes, hormones) which may detoxify contaminants (Lyons et al., 2015). This may lend some weight to an argument that an endogenous (indirect) interaction occurred in our experiment, in that copper was present at larger concentrations in mussels exposed to intermediate levels of chlorpyrifos than it was in all other stressor treatments ( Figure 2 ). Indeed, the potential for exogenous and endogenous interactions has been shown for other habitats and other stressors and those interactions also vary between being synergistic or antagonistic. An example of evidence of a synergistic exogenous interaction is one where the toxicity of certain polycyclic aromatic hydrocarbons is known to be enhanced on exposure to sunlight; in particular ultraviolet radiation (Oris and Giesy, 1985). Furthermore, the enhanced toxicity of polycyclic aromatic hydrocarbons has been shown to cause lethal effects to organisms at concentrations an order of magnitude smaller than those that are toxic without ultraviolet radiation (Ankley et al., 2010).

In contrast, examples of antagonistic exogenous interactions can occur with a range of metals in water where the toxicity and bioavailability of the metal are thought to depend on abiotic factors such as pH, temperature, salinity, “minerality-hardness” and concentrations of organic matter (Kiaune and Singhasemanon, 2011).

Such endogenous interactions can be even more complex because they may depend on biochemical and physiological processes in an organism (Lyons et al., 2015). For a range of aquatic species exposed to combinations of increased temperature and different chemical stressors there have been a mix of antagonistic or synergistic interactions. In particular, the direction of the interaction was dependent on not only the identity of the chemical stressor and the intensity of warming, but also the ability of either stressor type to rapidly activate a stress response; for example, the synthesis of “heat shock” proteins, which have also been shown to play a role in reducing toxicity of chemical stressors (reviewed in Holmstrup et al., 2010). Additionally, increased temperature combined with a non-ionic organic compound antagonistically effected freshwater sponges (Müller et al., 1995), but copper and warmer temperatures interacted synergistically in zebra mussels (Rao and Khan, 2000).

The interactive effect of copper and chlorpyrifos on the viability of hemocytes was not observed with filtration ( Figure 4A ). In fact, for filtering by assemblages our analysis indicates there were two non-additive interaction effects; a synergistic interaction as predicted under the additive model (Folt et al., 1999) and an indication that an interaction as predicted under the simple “comparative/dominance model” occurred (Folt et al., 1999; Côté et al., 2016). Firstly, there was a synergistic interaction at small to intermediate concentrations of chlorpyrifos and large concentrations of copper which reduced the filtering rate of assemblages compared to the controls that were not exposed to either stressor. Secondly, at larger concentrations of chlorpyrifos the addition of copper (across all concentrations) did little to change the effect of the chlorpyrifos, which in this instance determined the effect across the ‘response-surface’ ( Figure 4 ). Indeed, one possible explanation for the synergistic effect of stressors on the filtering by assemblages in our experiment, could be related to the synergistic increase in copper accumulating in mussels at these concentrations ( Figure 2 ). As explained above the mode of action for both these contaminants can be different, yet despite those differences the individual effects can combine to overwhelm the organisms at particular concentrations of both stressors and hence reveal synergistic effects or ‘surprises’ (Paine et al., 1998) at higher levels of biological organization. Previous experiments have shown that when mussels are exposed to large doses of chemical stressors, they can close their valves (Kramer et al., 1989), which may explain the increased copper in mussels exposed to small rather than larger concentrations of chlorpyrifos. Chlorpyrifos can be more toxic at smaller doses than copper (Serrano et al., 1997). Subsequently the effects of chlorpyrifos may reduce the organism’s ability to detect larger doses of copper and hence may ingest and store more copper when exposed to copper combined with small to intermediate doses of chlorpyrifos ( Figure 2 ), which is then translated to greater rates of filtering by assemblages.

Beds of bivalves are some of the most important and dominant habitat-forming species in aquatic systems, as they provide living space and favorable conditions for other species, which facilitates greater biodiversity (Norling and Kautsky, 2007). Mussel beds have diverse assemblages (O’Connor and Crowe, 2007), and their filtering can play a major role in assimilating and recycling sediment, nutrients and energy from the water-column (Newell, 2004). They also provide an important source of food for scavengers/predators (e.g. crabs, starfish, birds, humans) and underpin a vast array of both aquatic and terrestrial food webs (Paine, 1966; Norling and Kautsky, 2007). As such, the negative effects of stressors as described in the current experiment could have serious implications for the functions of mussels beds. In fact, loss of dominant species due to multiple stressors from higher trophic levels is predicted to have a greater impact than loss of a species from other trophic levels and may change the productivity and other useful functions of organisms and assemblages (Hughes and Connell, 1999; Coors and De Meester, 2008; Kendrick et al., 2019).

In this experiment we found possible evidence of a simple “comparative/dominance model” interaction between the stressors tested (Folt et al., 1999). In these types of interactions, the interaction is dominated by the effect of a single stressor (Folt et al., 1999). In these interactions it would appear that at sufficiently large doses/concentrations of one stressor it makes no difference to add the other stressor. Whilst this sort of outcome is valuable to determine managerial approaches for two chemical stressors, as given this result a manager or policymaker could try to remove the worst stressor from the ecosystem. If, however, the more harmful stressor is one that would be expected to occur under climatic change such as increased temperature, acidification, etc., it may not be possible to modify its effects (Harley et al., 2006; Russell et al., 2009; Cheng et al., 2015). These sorts of scenarios could be examined using a ‘response-surface’ design and provide useful information to managers and policymakers who wish to avoid harmful multiple stressor scenarios, something that has previously been called for in the literature (Kroeker et al., 2017; Varaksin et al., 2018), but remains rare (but see Rindi et al., 2022 for experimental approaches on biofilms) due to the complexity of the task and remains a knowledge-gap (Orr et al., 2020; Jackson et al., 2021).

Our results highlight that non-additive effects may be even more complex and may be even more common than previously thought (Crain et al., 2008; Darling and Côté, 2008; Holmstrup et al., 2010), particularly given that those interactions (or lack thereof) can vary substantially across different response variables and across varying concentrations of the stressors. In the current experiment, if we had used a binary or two factor by two level experimental design, we may have reached erroneous conclusions in terms of the existence, strength and direction of the interactions between stressors. While our results show the value of using ‘response-surface’ designs to examine the effects of multiple stressors across a range of intensities, these designs are limited by the number of stressors and organismal responses that can be investigated. So to ensure causal inferences gained from response-surface experimental designs are environmentally relevant will require that the stressors and responses (biological and ecological) chosen are based on those found to correlate in environmental surveys. Using linked surveys and experiments in this way is crucial to understand cause-effect relationships between stressors and biota is crucial to restoring, rehabilitating and preserving biota and their habitats (Underwood, 1980; Underwood and Peterson, 1988; Burton and Johnston, 2010).

The variation in the chemical, biological and ecological outcomes we observed with mussel beds exposed to different concentrations of copper and/or chlorpyrifos across levels of biological organization has important implications for monitoring and managing the impacts and risks of priority pollutants to the ecosystem. As such depending upon the assay used there is a possibility of concluding there was no effect, when in fact based on our findings there was a substantial effect. This potential for ambiguity across those levels of organization further highlights the need for more careful selection of assays and experiments to determine the nature, magnitude and mechanism of impacts of multiple contaminants and other stressors to the ecosystem (Underwood, 1980; Underwood and Peterson, 1988; Adams, 2005; Martínez-Crego et al., 2010; Browne et al., 2015). That said, this ambiguity in the response variables may be driven in part by experimental procedures used to measure respiration and filtration as we cannot be completely certain that the chambers used were sealed fully against the pontoon-surface. Equally that ambiguity may to some extent be explained by differences in the mechanism and timing of toxicity (i.e. adverse-outcome pathway; Kramer et al., 2011; Browne et al., 2015) for the various combinations of contaminants at different levels of biological organization.

There is a growing consensus that many, if not all terrestrial and aquatic habitats are impacted in some way by multiple anthropogenic stressors (Halpern et al., 2015; Lyons et al., 2015). As such there is a pressing need to understand and quantify those effects, if we are to preserve biodiversity and the valuable services they provide (Adger, 2006; Crain et al., 2008; Clements and Rohr, 2009). Experiments, such as the one described here, have the potential to unravel some of the observed complexities surrounding the effects and interactions of multiple stressors. Potential improvements would be to include more response-variables, particularly those at levels of biological organization that may be dampened by inter- and intra-specific interactions (Harley et al., 2006; Crain et al., 2008; Bennett et al., 2015). In addition, experiments incorporating larger temporal and spatial scales could test if the types of interaction vary or are general through time and/or space, or indeed may vary due to the timing/frequency and duration of exposure to the stressors themselves. Understanding how multiple stressors interact to impact the ecosystem remains a challenge but this challenge needs to be met to identify and manage harmful anthropogenic activities so that they have smaller ecological impacts.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

PB, MB, LB-C and TC conceived the study. PB and MB set up and ran the experiment, while DL assisted PB, MB and others with field-measurements on the final day of the experiment. PB organized the database. PB, LB-C and DL performed the statistical analysis, PB and DL created the figures. TC supervised PB’s PhD of which this study was part of. PB wrote the first draft of the manuscript. All others contributed to manuscript revision and approved the submitted version.