Copper mediates life history responses of Daphnia pulex to predation threat

Introduction

Species within ecosystems are under constant pressure from a mixture of biotic and abiotic challenges that represent a complex mosaic of selective pressures. This is particularly true for freshwater systems where the majority have been modified and continue to be influenced by human activity including agriculture, sewage, and mining (Dudgeon et al., 2006; Dudgeon, 2019; Reid et al., 2019). As a consequence, biodiversity in freshwater habitats is at high risk: while such habitats constitute less than 1% of the earth’s surface, research estimates that freshwaters harbor approximately 6% of all species (Dudgeon et al., 2006).

Anthropogenic modifications to these ecosystems occur in parallel to a mosaic of species interactions like predation and disease that determine the traits, distribution and abundance of species; anthropogenic stressors are likely to interact with these interactions (Relyea, 2005; Moe et al., 2013; Katzenberger et al., 2014; Reid et al., 2019). Multi-stressor research over the past decades has extensively broadened our understanding of how different stressors combine to impact on species (Folt et al., 1999; Holmstrup et al., 2010; Altshuler et al., 2011; Jansen et al., 2011) and highlights physiological and life-history traits as “targets” of stressors that give rise to modified patterns of population growth, abundance and distribution. This research provides a solid foundation to address the question of when and in what traits might we detect synergistic or antagonistic interactions between anthropogenic effects and the impacts of natural “threats” like predation and disease.

Here we focus on potential interactions between heavy metals and predation risk detected in physiology, morphology and life history traits of a model organism in both ecology and ecotoxicology, Daphnia pulex. The nature of exposure to heavy metals in freshwater habitats is now largely sub-lethal. This exposure paradigm aligns with the extensive history of research on the sub-lethal effects of predation as a natural form of stress: sub-lethal effects are typically studied in the framework of phenotypic plasticity where organisms’ physiology, morphology and life history take on values that depend on the environmental conditions (Ghalambor and Martin, 2002; Auld and Relyea, 2011; Hawlena et al., 2011a,b).

Plastic response to sub-lethal predation risk are typically defined as adaptive (Ghalambor and Martin, 2002; Auld et al., 2010): changes in physiology, morphology and life history are aligned with increasing survival and maintaining reproduction in the face of mortality risk at specific ages. In contrast, sub-lethal responses to heavy metals are often still defined as acute or chronic negative effects reflecting a history of research on toxicity (De Schamphelaere et al., 2005; Moe et al., 2013; Brockmeier et al., 2017; Stubblefield et al., 2020). Yet, the response to, for example, starvation by digestive inhibition by a metal, may elicit classic, adaptive life history responses to reduced energy intake (Ferrando and Andreu, 1993; De Schamphelaere et al., 2007). Thus, in order to understand how natural and anthropogenic stressors interact, it is vital to integrate knowledge about how physiology, morphology and life history respond, via plasticity, to multiple types of sub-lethal stress targeting these traits.

Predation risk versus copper exposure: Sub-lethal responses

One of the most ubiquitous features of natural communities is the risk of predation. Research centred on the waterflea D. pulex is one of the best studied and popular examples of predator induced phenotypic plasticity (Tollrian, 1993; Tollrian, 1995; Beckerman et al., 2010; Dennis et al., 2011; Beckerman et al., 2013; Carter et al., 2017). Decades of D. pulex research has documented an astonishing repertoire of inducible defenses reflected in predator cue induced changes to morphology and life history (Agrawal, 1999; Riessen, 1999a; Beckerman et al., 2010; Dennis et al., 2011; Riessen, 2012; Carter et al., 2017): Exposing D. pulex to predator cues induces both morphological defenses and substantial shifts in the size and age at maturity and the allocation of resources to growth vs. reproduction. There is immense genetic variation and among population variation in induced defenses and evidence for genotypes with little or no response (low-responders) and for genotypes with very strong responses (high-responders; Lind et al., 2015; Spitze, 1992; Weber and Declerck, 1997). Variation in life history arises from variability in size selective predation where small size selective predation favors somatic growth over reproduction in the prey, reflected in a later reproduction at a larger size, while large size selective predation favors reproduction over growth and leads to an early maturation at smaller sizes (Riessen, 1999a; Beckerman et al., 2010; Riessen and Gilbert, 2019).

Heavy metals remain one of the most prevalent anthropogenic pollutants in the aquatic environment (e.g., Häder et al., 2020), with copper (Cu) being considered one of the most toxic metals for aquatic ecosystems (Nor, 1987). While Cu is an essential resource for organisms, particularly in the context of redox reactions, elevated Cu concentrations in freshwater habitats due to agricultural, municipal, and industrial discharges remain of concern world-wide (e.g., Rader et al., 2019). While concentrations in the wild are often not toxic, pervasive exposure to sub-lethal concentrations of copper continues. These sub-lethal concentrations strongly affect foraging and the assimilation of nutritional resources, leading to reduced energy intake and increased metabolic costs (De Coen and Janssen, 1997; Wilding and Maltby, 2006). These alterations to foraging typically translate into reduced reproduction and somatic growth associated with smaller size and later ages at maturity, patterns that align with life history plasticity theory about the impacts of starvation (Van Noordwijk and de Jong, 1986; Boggs, 1992).

The work outlined above highlights the central role of resource allocation to growth and reproduction, leading to distinct life history outcomes in size and age at maturity and reproduction, for both predation and Cu. Despite extensive research on each, studies on their interaction are scarce and potentially limited by the number of genotypes and traits considered (e.g., Clements et al., 1989; Hunter and Pyle, 2004; Mirza and Pyle, 2009; DeMille et al., 2016). The life history theory details above, however, provide a framework for predicting how these two stressors may interact. We summarize these predictions in Figure 1, providing the basis for a formal test for synergistic or antagonistic effects of the two stressors. The key information in this table is that life history theory suggests the potential for some synergistic (age at maturity) and some antagonistic (size at maturity, somatic growth rate) effects of multiple stressors.

Figure 1. Candidate predictions for plastic responses to copper and predation risk by two classes of D. pulex genotypes: class one has a strong life history response, but a weak morphological response to predation risk (“low responder”) while class two is the opposite (“high responder”). As discussed in the text, copper is expected to reduce size, increase age and decrease growth rate due to reduced energy intake arising from disrupted digestive functioning (acquisition and assimilation). In contrast, predation risk is expected to increase size and decrease age at maturity and increase growth rate. This reveals the potential for some synergistic (age) and some antagonistic (size, growth rate) effects of multiple stressors.

To address the predictions in Figure 1, we exposed replicates of six different genotypes of D. pulex to Cu and predation risk in a response surface statistical design comprising four concentrations of each stressor (see methods). We focus our assessment on five key traits (see Figure 1): induced morphological defenses, size at maturity, age at maturity, somatic growth rate and survival until maturity. We explore the potential for these interactions between predation risk and Cu in six genotypes, three with a known “strong” morphological versus life history response to predation risk and three with a known “weak” morphological versus life history response. This design provides a robust framework, grounded in life history theory, to differentiate between additive and non-additive effects of the stressors among five traits. The data also allow us to propose a hypothesis that there are different strategies involving multiple traits when dealing with multiple stressors.

Materials and methods

Daphnia gentoypes

Six D. pulex genotypes (D8.6A, D8.7A, Cyril, Chard, C14, LD3.3) were collected from different field populations (D = Dorset, C = Crabtree, LD = Lake District, United Kingdom; for details on clones and sites see Supplementary Table S1) in 2010 and maintained as iso-female lineages in artificial hard water (ASTM; (ASTM, 1988)) with seaweed extract (marinure; Wilfrid Smith Std., Northans, United Kingdom; 0.036 μl ml⁻¹) under standard conditions: 15 animals L⁻¹ ASTM were fed three times a week with Chlorella vulgaris (2 × 10⁵ cells ml⁻¹; >1.5 mg carbon L⁻¹) and maintained under 16:8 h light:dark conditions at 20 ± 1°C. Genotypes were molecularly defined by a set of 11 microsatellite markers (Reger et al., 2018).

Multi-stressor environment

We generated predator cues from frozen midge larvae following established protocols (Tollrian, 1995; Beckerman et al., 2010). Briefly, homogenized midge larvae extracts were filtered, followed by solid-phase extraction using a C₁₈ column (Agilent) to recover the active compounds that generate strong morphological and life history responses in D. pulex.

Copper (Cu) stock solutions were prepared by dissolving copper (II) chloride dihydrate salt (Fisher Scientific, United Kingdom) in distilled water, which were diluted with ASTM medium to create experimental exposure conditions. Realized Cu concentrations were tested by collecting medium samples from experimental exposure treatments throughout the study and using ICP-MS; realized concentrations did not deviate from nominal concentrations (Cu: R² = 0.99, F = 1913, p < 0.002; see Sadeq and Beckerman, 2019).

In order to produce a realistic fine-scale gradient of both stressors (i.e., risk of predation and copper exposure), we exposed animals to each of 16 test conditions, arising from four concentrations of each stressor (predation cue: 0, 0.1, 0.25 and 0.5 μl mL⁻¹; copper: 0, 5, 10, 25 μg L⁻¹) and all combinations. All used concentrations are environmentally relevant and were chosen based on previous research (see, e.g., Lind et al., 2015; Beckerman et al., 2010; Dennis et al., 2011; Sadeq and Beckerman, 2019, 2020).

Experimental exposure

We performed a classic life table experiment with n > 3 replicates of each genotype in each treatment condition to collect data on five traits. Prior to experiments, animals were kept under standard conditions (food = 2×10⁵ cells ml⁻¹ Chlorella vulgaris; light:dark = 16:8 h, no predation risk) for three generations but fed daily. For experiments, two to six gravid Daphnia, carrying embryos in E₃ stage (~18 h before parturition; sensu Naraki et al., 2013), were placed in individual jars containing 50 ml ASTM, algae (2 × 10⁵ cells ml⁻¹ Chlorella vulgaris), liquid seaweed extract (marinure; Wilfrid Smith Std., Northans, United Kingdom; 0.018 μl ml⁻¹) and a combination of Chaoborus predator cue and copper. After parturition, two neonates were randomly selected from each mother per treatment (n = 4 biological replicates in total) and placed individually in 50 ml glass vials containing the same medium as their maternal environment. Until maturity, animals were photographed daily (MZ6 modular stereomicroscope, Leica; EOS 350D camera, Cannon) and subsequently transferred to a new glass vial containing fresh media, algae (2 × 10⁵ cells ml⁻¹ Chlorella vulgaris) and respective environmental stressor cues. To assess predator induced defense morphologies from images, appropriate magnification was used. Due to the large amount of experimental exposures, phenotypic assessments were carried out in a split-block experimental design across 4 batches with all six clones phenotyped concurrently and experimental exposure jars randomized across the rearing set-up and evenly split between treatment groups.

Data collection

Using daily photographs, we captured data on five response variables: induction of morphological defenses at instars 1–3, size at maturity, age at maturity, somatic growth rate, and survival rate until maturity. Morphological induction, based on the presence of a pedestal and spikes on the neck of the Daphnia, was manually scored using a previously defined scoring system (Tollrian, 1995; Hammill and Beckerman, 2010): pedestals were classified as absent (score = 0), small (score = 30) and large (score = 50). Spikes in the nuchal area, later referred to as neckteeth, were counted individually and provide an additional score of 10 each. For further analyses, we used the maximal induction that occurred over the course of the first three instars. Age and size at maturity were estimated as size and age of adult animals on the day neonates first appeared in the brood pouch (max 13 days). For animal size, we manually measured animal length from the tip of the head to the base of the tail spine using the image analysis software ImageJ (Schneider et al., 2012). Somatic growth rate was calculated as (size_Maturity − size_Initial)/age_Maturity.

Detecting non-additivity of predation risk and copper exposure

To assess the contribution of genotype, treatment, and their interaction on phenotypic traits, we fit a classic response surface model within a mixed effects framework using the lme4 package for R (Bates et al., 2015). The full model to detect non-additive effects and characterize the shape of the response surface takes the form of response ~ PredationRisk + PredationRisk² + Cu + Cu² + PredationRisk*Cu + (1|Genotype) where (1|Genotype) represents the random effect of genotype. The fixed effects, response surface component is a standard statistical model with applications in ecology and ecotoxicology that can detect planar and several non-linear “response” surfaces (e.g., peaks, troughs and valleys) driven by two covariates (Jonker et al., 2005; Scholes et al., 2005; Sadeq and Beckerman, 2020). The covariates are modeled as polynomial terms such that the effects of predation risk and Cu can be linear or non-linear and this is combined with an explicit interaction term, which together with the non-linear terms, allows for various shapes of an interaction surface to emerge if justified by the data. For survival data, we fit the same model in a cox proportional hazard framework with a frailty term for clone. Survival data are presented as risk scores that quantify the risk of death: higher numbers (scores) indicate a higher risk of death and thus a reduced probability of survival. To evaluate the response surface model and evidence for non-linear responses and interactions for each trait (see Supplementary Table S2) we employed Type II sums of squares implemented by the Anova() function in the car package for R (Fox and Weisberg, 2019).

Strategies for managing two stressors: A hypothesis

To develop the hypothesis that clones with strong versus weak morphological defenses (clone type) to predation risk respond fundamentally differently to combined exposure to predation and copper, we first coded the data to ask whether a predation x copper interaction varied by clone type. To estimate this three-way interaction, we fit linear models for each trait: response ~ PredationRisk + PredationRisk² + Cu + Cu² + PredationRisk * Cu * CloneType. This model differs from the model above by replacing the random effect of clone with a fixed effect of clone type that divides the clones into two groups: weak responders and strong responders (re the development of defense morphologies).

The presence of a three-way interaction (i.e., predator cue * Cu * clone type), for life history traits in particular, would suggest that intrinsic levels of morphological defense to predators might drive trade-offs between life history and morphological responses that lead to conditional outcomes of exposure to predation risk and Cu. For survival data, we evaluated the same three-way interaction via a cox proportional hazard model.

To visualize and further evaluate these potential trade-offs that might support a hypothesis that there are different strategies for dealing with multiple stressors, we plotted the relationships between three life history traits and induced morphology, at the extreme predation and copper treatments, i.e., high copper concentration & no predation cues (0 μl mL⁻¹ predation cue; 25 μg L⁻¹ copper), high predation cue concentration & no copper (0.5 μl mL⁻¹ predation cue; 0 μg L⁻¹ copper), and high predation cue concentration & high copper concentrations (0.5 μl mL⁻¹ predation cue; 25 μg L⁻¹ copper) for the strong and weak responders in contrast to control conditions (i.e., no predator cue and no copper). Specifically, this involves plotting the bi-variate mean and SE of trait combinations (i.e., life history and morphology) for each group and each treatment. Using this visualization, we can characterize on three patterns. First, are clones that possess weak morphological response responding weakly in bivariate trait-space compared to clones that show strong morphological response? This would be evidenced by smaller vectors between control and non-control treatments. Second, do clones that possess weak morphological response react to copper and predation (single threats) differently? This would be evidenced by the orientation of the vectors between control and each stressor. Finally, do clones that possess weak morphological response react differently to the combined threats? This would be evidenced by different lengths and orientation of the vector between control and dual-threat conditions.

Finally, to add some quantitative support for the hypothesis that there are different strategies, we apply a Procrustes analysis (procOPA function from the shapes package in R; see Dryden, 2021) to the convex hull that forms from connecting the response of multiple traits to the multiple stressors. This analysis provides estimates of the change in location, size and shape of the hulls in each panel of the visualization, along with an estimate of root-mean-squared deviation.

All analyses were made using R 4.0 (R Core Team, 2020). Visualization and data manipulation was supported by the following packages: ggplot2 (Wickham et al., 2016), cowplot (Wilke, 2020), data.table (Srinivasan, 2021), and viridis (Garnier et al., 2018).

Results

Overall we found that, in all life history traits, and for each stressor, observed responses typically aligned with classic life history theory about predation risk and heavy metal induced starvation risk. However, we found no evidence that predator induced defense morphologies varied by copper. We found evidence in several life history traits that the effect of predation risk was modified by the heavy metal. We found evidence only in age at maturity that the interaction between predation risk and copper varied with the designation of clones as high vs. low morphological responders. Despite this, we see evidence in bi-variate trait space that these intrinsic phenotypic details might underpin different strategies. Here we report on each analysed trait (Figures 2–6; Supplementary Figures S1–S5) with a focus on whether and how the effect of predation risk on traits is modified by Cu concentrations. We then present insight into potential trade-offs associated with intrinsic levels of induced defenses (Figure 7; Supplementary Figure S6).

Figure 2. Morphological induction in two classes of D. pulex genotypes. (A) Contour plots of the average fixed effects portion of the model showing the overall morphological changes to copper and predation in “high inducers” and “low inducers.” (B,C) Average effects of risk of predation (B) and sub-lethal copper concentrations (C) on morphological induction when copper and risk of predation, respectively, are absent (e.g., references for each treatment). X-axes present predator cue (A,B) or copper concentrations (C), and y-axis denote copper concentrations (A) or morphological induction estimates (B,C).

Figure 3. Size at maturity in two classes of D. pulex genotypes. (A) Contour plots of the average fixed effects portion of the model showing the overall response (here: size at maturity) to copper and predation in “high inducers” and “low inducers.” (B,C) Average effects of risk of predation (B) and sub-lethal copper concentrations (C) on size at maturity when copper and risk of predation, respectively, are absent (e.g., references for each treatment). X-axes present predator cue (A,B) or copper concentrations (C), and y-axis denote copper concentrations (A) or size at maturity estimates (B,C).

Figure 4. Somatic growth rate in two classes of D. pulex genotypes. (A) Contour plots of the average fixed effects portion of the model showing the overall response (here: somatic growth rate) to copper and predation in “high inducers” and “low inducers.” (B,C) Average effects of risk of predation (B) and sub-lethal copper concentrations (C) on somatic growth rate when copper and risk of predation, respectively, are absent (e.g., references for each treatment). X-axes present predator cue (A,B) or copper concentrations (C), and y-axis denote copper concentrations (A) or somatic growth rate estimates (B,C).

Figure 5. Age at maturity in two classes of D. pulex genotypes. (A) Contour plots of the average fixed effects portion of the model showing the overall response (here: age at maturity) to copper and predation in “high inducers” and “low inducers.” (B,C) Average effects of risk of predation (B) and sub-lethal copper concentrations (C) on age at maturity when copper and risk of predation, respectively, are absent (e.g., references for each treatment). X-axes present predator cue (A,B) or copper concentrations (C), and y-axis denote copper concentrations (A) or age at maturity estimates (B,C).

Figure 6. Risk of death in two classes of D. pulex genotypes. (A) Contour plots of the average fixed effects portion of the model showing the overall response (here: risk of death before maturity) to copper and predation in “high inducers” and “low inducers.” (B,C) Average effects of risk of predation (B) and sub-lethal copper concentrations (C) on risk of death when copper and risk of predation, respectively, are absent (e.g., references for each treatment). X-axes present predator cue (A,B) or copper concentrations (C), and y-axis denote copper concentrations (A) or risk of death estimates as risk scores (B,C).

Figure 7. Clone type specific responses in key life history traits. Plots contrast mean morphological induction vs. observed responses in key life history traits: size at maturity (A), somatic growth rate (B), and age at maturity (C). Black points: low responders; red points: high responders. Error bars indicate standard errors of mean values. Note: graphs show responses at extreme points in the response surface, i.e., control conditions (i.e., no predator cue and no copper) (a), high predation cue concentration & no copper (b), high copper concentration & no predation cues (c), and high predation cue concentration & high copper concentrations (d).

Induction of morphological defenses

We found no evidence that the effect of predation cues on defense morphologies varied by Cu (Figure 2A; Supplementary Figures S1B,C; Supplementary Table S2). Copper had no effect on induced defenses (Figure 2C; Supplementary Figure S1E; Supplementary Table S2, main effect “Cu”). We observe levels of morphological defense predicted from theory and previous work (Figures 1, 2A,B; Supplementary Figure S1A–D; Supplementary Table S2): Three clones (i.e., C14, Chard and LD33) developed pronounced defense morphologies in the nuchal area, with levels of induction increasing with increasing risk cue concentrations while the three other clones (i.e., D8.6A, D8.7A and Cyril) showed low levels of morphological induction that were independent of the tested cue concentration. In these “low responders,” morphological responses rarely involved the development of a distinct pedestal but the development of a single to a few small neckteeth. Hereafter, we will refer to these two groups of clones as “high-responders” and “low-responders,” respectively.

Size at maturity

Our data indicate that copper mediated the impacts of predation risk on size at maturity (Supplementary Figure S2; Supplementary Table S2). As copper concentrations increase, the increase in size at maturity caused by predation risk is “flattened,” largely because Cu reduces size at maturity more at low predation risk than at high (Supplementary Figure S2B). This finding suggests that the simultaneous exposure to both stressors induces a non-additive, “antagonistic” effect as predicted in Figure 1. In line with our predictions (see Figure 1), the risk of predation increased size at maturity and copper reduced size at maturity (Figures 3B,C; Supplementary Figure S2D,E; Supplementary Table S2).

Somatic growth rate

Our data reveal that the effect of predation risk on somatic growth rate is mediated by Cu (Figure 4A; Supplementary Figure 3B,C; Supplementary Table S2): Without Cu exposure, growth rate decreased with increasing predator cues. This pattern is essentially reversed at high levels of Cu where growth is strongly reduced at low predation risk (Figure 4A; Supplementary Figure S3B). This aligns with prediction in Figure 1 for an antagonistic effect. In line with predictions (see Figure 1), predation risk had negligible impacts on somatic growth rate while copper caused a strong reduction (Figures 4B,C; Supplementary Figure S3D,E; Supplementary Table S2).

Age at maturity

Our data show that at low concentrations of Cu, there is a slight increase in age at maturity with increasing predation risk. However, at high Cu, this is reversed because Cu increases the age at maturity at low predation risk more than it does at high predation risk, where age is already higher (Figure 5A; Supplementary Figure S4B; Supplementary Table S2). This finding suggests that the simultaneous exposure to both copper and predator cues is synergistic, in line with predictions from Figure 1. Partially in line with our predictions (see Figure 1), we found that predation risk generated negligible effects on age at maturity while copper dramatically extended the age at maturity (Figures 5B,C; Supplementary Figure S4D,E; Supplementary Table S2).

Survival rates until maturity

Our data indicate that the effect of predation risk on survival until maturity is mediated by Cu (Figure 6; Supplementary Figure S5; Supplementary Table S2). An increased risk of death (here: increase in risk scores - darker to lighter colors) caused by predator cues at low levels of Cu is reversed at high levels of Cu (Supplementary Figure S5A; Supplementary Table S2) as high Cu concentrations cause a strong reduction on survival at low predation risk. This aligns with predictions in Figure 1 for an antagonistic effect.

Copper concentrations used in this study were chosen based on survival estimates after exposing Daphnia for 24 h and 48 h to diverse Cu concentrations ranging from 5 μg L⁻¹to 450 μg L⁻¹ (Sadeq and Beckerman, 2019). While these data indicated that there was no difference in the response of neonates from all tested clones to an acute exposure to diverse sub-lethal copper concentrations, in our experiments survival during development was strongly affected by exposure to the highest chosen copper concentration (Figures 6A,C; Supplementary Figure S5A,D; Supplementary Table S2). Predation risk, however, had only minor effects on survival (Figures 6A,B; Supplementary Figure S5A,C; Supplementary Table S2). Both findings support our above stated expectations (see Figure 1).

Strategies for managing two stressors: A hypothesis

We only found evidence for the proposed three-way interaction among copper, predation and clone type in age at maturity. For all other traits, the effects of predation and/or copper, but not their interaction, varied by clone type (e.g., two-way interactions; see Supplementary Table S2).

To further explore whether there might be different strategies in managing threats by low versus high responders to predation risk, we show in Figure 7 the bivariate relationship between three life history traits and the morphological response (bivariate means ± SE) at the extreme predation and copper treatment levels compared to control conditions, colored by clone response group.

We posed three questions in the methods and we see evidence for all three: low responders respond less in multiple traits to the single and combined threats; low responders have a unique response to each stressor, partially driven, we hypothesize, by the hard boundary of zero morphological response (the predation vector cannot go below 0); and low responders have a vector of response to combined stressors that is more biased toward responding to predation, even though their levels of response to threats are very small. These results combine to reflect clear difference in size and location of the response “hull” and subtle difference in the “shape” of the hulls.

Finally, to further support the development of the hypothesis, we applied a Procrustes analysis to compare the two convex hulls in each panel of Figure 7 (see Supplementary Figure S6). This resulted in estimates of the size and shape of the hulls in each panel along with an estimate of root-mean-squared (rmsq) deviation. The location (i.e., the average morphology-life history relationship) varies dramatically (Figure 7; Supplementary Figure S6). The strong inducers have a much larger total response (all Δ centroid-size >32). The scaling of all high to low responders had large values (all >3) and the rmsd for all comparisons is large (>13). Together Figure 7 and the shape analysis (Supplementary Figure S6) help define the hypothesis that low and high responders operate in the morphological-life history trait response space differently.

Discussion

Natural freshwater habitats are home to myriads of species, harboring a large fraction of Earth’s biodiversity, while at the same time providing vital ecosystem services including drinking water, food and recreation (Dudgeon et al., 2006; Dudgeon, 2019). However, organisms living in such habitats experience a mosaic of natural threats (e.g., predation and disease) who’s impacts on life history are potentially modified by anthropogenic stress (e.g., Hendry et al., 2017). The ability to detect and adequately respond to this array of environmental cues is a vital feature of organism biology that directly affects survival and fitness, and determines their abundance and distribution. Keystone species’ responses are of particular interest because such interactions can trigger substantial cascading consequences for an entire community or ecosystem (Mills et al., 1993).

In particular, metal pollution of freshwater habitats has drawn much attention due to their persistence and deleterious effects on aquatic food chains. Even though there has been considerable exploration into the toxicity of different heavy metals, including copper, and their effect on aquatic organisms, research is typically focused on single stressors, single genotypes or single concentrations that are beneficial to a testing environment but not reflective of natural communities (e.g., Adams et al., 2000; LaLone et al., 2017). Such research thus not reflects “real world scenarios” despite the urgent need to do so (Dudgeon et al., 2006; Dudgeon, 2019; Reid et al., 2019).

Here, we evaluate how a keystone species of freshwater communities world-wide, D. pulex, responds to a mixture of the natural risk of predation when exposed simultaneously to an anthropogenic stressor - the heavy metal copper. In order to understand how anthropogenic impacts interfere with or mediate natural environmental cues, we specifically asked whether both stressors combine to affect key morphology and life history traits in an additive or interactive manner. Our analyses are framed around ecological and ecotoxicological theory associated with resource limitation and similar theory about predation risk by small size selective predators experienced by D. pulex (see Figure 1). Our data are unique in assessing the impacts of gradients of the two stressors on five traits and evaluating responses among six genotypes.

Clonal variation in anti-predator responses

Supporting previous findings (Dennis et al., 2011; Carter et al., 2017), we found morphological response to predation risk exposure to vary among genetically distinct clones: three clones responded to predator cues with strong morphological induction, while the other three clones showed minor morphological changes, if any at all. This finding parallels well established evidence of genetic variation in induced defenses (Lind et al., 2015; Spitze, 1992; Weber and Declerck, 1997) and that differences in morphological defenses reflect natural variability that is likely driven by local adaptation and potentially different strategies for managing predation risk (e.g., Reger et al., 2018). Despite significant differences in morphology, all clones showed strong predator cue driven responses in life history that are well established theoretically and empirically for daphnids (Tollrian, 1993; Tollrian, 1995; Hammill et al., 2008; Beckerman et al., 2010; Hammill and Beckerman, 2010): we found animals’ size at maturity to increase upon exposure to predator cues. However, this size increase was not due to increased somatic growth rates, but resulted from morphologically high inducing clones being larger in general. Contrary to previous reports (Riessen, 1999b; Beckerman et al., 2010), we also found no significant effect of predator cues on the age at maturity for the tested clones.

Copper mediates life history responses but not morphology

Our results across multiple clones reveal no evidence that the effect of predation risk on morphology was mediated by copper: clones that strongly responded to predator cues show the same extent of morphological induction when exposed to copper as well. The ability to maintain the development of anti-predator responses in the face of copper suggests that it does not directly interfere with the signal transduction pathways underpinning the induction of anti-predator defense. Our finding partly contradicts previous reports (Hunter and Pyle, 2004; Pyle and Mirza, 2007; DeMille et al., 2016) that found a strong interaction effects between the heavy metal and predator cue induced morphological changes. However, some of these analyses (Hunter and Pyle, 2004; Pyle and Mirza, 2007) reported on the responses of single genotypes and may have focused - based on the reported animal sizes - on the morphological inductions at 1^st instars and not 2^nd instars, or the maximal induction over the first few instars, or solely the number of neckteeth (ignoring the presence of a pedestal). This work may have thus not assessed the full extent of morphological induction at the stages where it generates the most benefit. Further work from the same group (DeMille et al., 2016) later analysed potential effects of a combined exposure to the heavy metal and predation risk on a larger number of genotypes, but that study only reported on life history changes and not on morphological inductions, and thus does not allow for an adequate comparison.

We found copper to have a profound effect on life history traits, in line with previous findings (De Schamphelaere et al., 2007; Agra et al., 2011; Sadeq and Beckerman, 2019; Sadeq and Beckerman, 2020): exposure of animals to a gradient of copper concentrations alone resulted in a strong reduction of size at maturity and somatic growth rate, in addition to substantial increases in the age at maturity. Our data further indicate that an exposure to high levels of copper reduced the chances of survival until maturity. For our study, we can exclude the possibility that the observed responses to the heavy metal were due to clones being tolerant to the tested copper concentrations. In fact, while none of the here used copper concentrations induced any nominal mortality rates in neonates across 24- and 48-h test trials, exposure to the highest concentration of the heavy metal (i.e., 25 μg L⁻¹) during the full development from neonate to adult stages resulted in an increased risk of death. This result demonstrates that while metal concentrations in the environment may not be high enough to cause immediate toxic effects, sub-lethal concentrations can induce pronounced shifts in life-history that make species more susceptible to environmental threats, leading to a reduced overall fitness. A prolonged exposure to copper during the development likely results in an accumulation of internal copper concentrations, which may both impair cellular processes and require energy resources to be diverted toward distinct cellular detoxification processes (Bossuyt and Janssen, 2004). Our findings may also suggest that the digestive inhibition mode-of-action of copper induces starvation over the life cycle, with the reduced energy uptake responsible for the observed mortality in some of the genotypes (Ferrando and Andreu, 1993; De Schamphelaere et al., 2007).

Our analyses provide further evidence that copper modified the effect of predation risk on several key life-history traits. Overall, the observed responses on life history traits were in line with expected changes in life history by each stressor (see Figure 1) and allow differentiating between additive and non-additive effects of stressors.

Specifically, we found evidence of antagonistic effects for size at maturity, somatic growth rate and survival until maturity, when animals were exposed to both stressors simultaneously. A possible explanation for this, in line with theory and evidence outlined via Figure 1, is that two different allocation mechanisms - one tied to starvation and another to size selective predation - are offsetting each other in a non-additive manner (see Noonburg and Nisbet, 2005; Beckerman et al., 2010). In contrast, but expected, for age at maturity, we observed a non-additive, synergistic behavior of the two stressors, likely for the same reason as the allocation strategies are aligned here (see Figure 1).

Are there different strategies defined by low and high responders?

Based on theory that resource allocation plays a vital part in morphological and life history responses to stressors, we presented the hypothesis that high inducing Daphnia clones, that may allocate more resources to the defense morphology than to life history, might differ in their response to Cu and predation risk when compared to clones that show minor anti-predator defense morphologies. We note that our data are sufficient only to propose, not test this hypothesis.

Our response surface analyses, aiming at understanding the relationship between copper and predation over gradients of the two stressors, provide such evidence only for one of the tested life history traits: age at maturity. This result points to the potential role of age at maturity as a pivot point for managing multiple stressors and that allocation to morphology may trade-off with this life history trait that is key on it’s own and also part of somatic growth rate plasticity.

Despite this lack of a three-way interaction in many traits, our assessment of bivariate relationships between morphological defense and life history traits (Figure 7) - e.g., the potential trade-off between these traits - suggest that there are some substantive differences between clones that respond weakly or strongly to predation risk with morphology. For all three life history traits in Figure 7, and with evidence from a Procrustes analysis, the location, size and shape of the hull connecting control and treatment responses are different (see Supplementary Figure S6). We refrain from making strong inference here as there are only three genotypes in each category. Our data suggest that future efforts to characterize multivariate responses to multiple stressors might reveal strong patterns that define trade-offs among traits and threats. Such data may provide insight into whether strong and weak morphological defenses lead to altered allocation of the required resources for other predation risk responses and other cellular processes, including metal detoxification.

Conclusion

Overall, we found no evidence that the heavy metal copper mediates the effects of predator cue induced morphological responses in the tested D. pulex clones. However, our data indicate that copper mediates several key life-history responses to predation risk, with all observed responses supporting predictions from life-history theory. We also found evidence that univariate and multivariate analyses with high replication are likely necessary to effectively reveal how resource allocation to defenses against multiple stressors may trade-off via different strategies. Our study also emphasizes that if we are to move toward generalizing and predicting the effect of multiple stressors, experiments need to account not only for the diversity of traits that can be evaluated but also requires sophisticated approaches that accurately accommodate for key features of natural communities: (1) stressors do not occur alone and vary in concentrations; (2) there is genetic variation in multiple stressor responses; and (3) this genetic variation likely underpins local adaptation. These three features of real-world communities present challenges of predicting the effects of multiple stressors, but make imperative the need to do so. Here, we provide evidence that “the use of response surface models is an effective and powerful way to assess additive and non-additive effects, providing statistically based conclusions about combined dose responses” (see also Jonker et al., 2005).

Data availability statement

All scripts and code used for data analysis and plotting are available at https://github.com/beckerdoerthe/MultiStressorPlasticity.

Author contributions

DB: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, resources, visualization, writing - original draft, writing - review and editing. AB: conceptualization, funding acquisition, methodology, resources, writing - review and editing. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the European Union’s Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie grant agreement no. 841419 and by awards #BE 5288/2-1 and BE 5288/3-1 from the German Research Foundation. The authors also acknowledge the University of Sheffield Institutional Open Access Fund.

Acknowledgments

The authors would like to thank Allison Blake and Lynsey Gregory for technical support. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.983923/full#supplementary-material

References