From organismal physiology to ecological processes: Effects of nutrient enrichment and salinity variation in a freshwater ecosystem

Environmental stressors have profound implications for species, communities, and ecosystems by altering fundamental processes. With increasing human impacts on aquatic ecosystems, two main scenarios have been reported: (1) the spatiotemporal superposition of multiple stressors, leading to interactions among them and (2) intensifying environmental gradients, leading to threshold responses. However, studies designed to assess the effects of multiple stressors (e.g., pulse stressors) along environmental gradients (e.g., press stressors) are uncommon, and interactions between pulse and press stressors may cause abrupt changes in biological responses. We conducted a laboratory experiment to investigate the effects of osmotic stress along a nutrient enrichment gradient on a freshwater community composed of periphyton, microorganisms, and zebra mussels (Dreissena polymorpha). Our objectives were to (1) quantify the individual and combined effects of stressors, (2) delineate thresholds along the nutrient gradient (press) in the absence and presence of osmotic stress (pulse), and (3) test for interactions between the two stressors. We evaluated effects on metabolic rates in D. polymorpha and on microbial activity, as well as phototrophic periphyton biomass and physiological status. We observed interactions between the two stressors for metabolic rates in D. polymorpha and periphytic phaeopigments. In contrast, we found an individual effect of osmotic stress on microbial activity and chlorophyll a content. Thresholds were only identified in the presence of osmotic stress for metabolic rates in D. polymorpha. Our work highlights the importance of combining multiple stressors with environmental gradients and the need to consider multiple biological compartments when evaluating the impacts of stressors on ecosystems.

Increases in human population and activity are leading to sudden changes in natural habitats, populations, and ecosystems (Halpern et al. 2019). The pressures deriving from these activities, known as environmental stressors or drivers, can be defined as factors of natural or anthropogenic origin perturbing a biological component beyond its natural limit of variation and tolerance (Vinebrooke et al. 2004;Carrier-Belleau et al. 2021a). The origin of these stressors occurs at a range of spatial scales, from local anthropogenic activities (e.g., land use, coastal development) to anthropogenically induced global climate change (e.g., warming, acidification; Hewitt et al. 2016). Stressors may impact different levels of biological organization (Gissi et al. 2021) and ultimately affect ecosystem functioning and services (Balmford et al. 2011).
Due to the intensified human footprint in natural ecosystems, the accentuation of environmental gradients has been observed spatially and temporally (Scheffer et al. 2001). While individuals and ecosystems may respond linearly to a gradient, they can also reach thresholds that correspond to a level of environmental pressure beyond which there is a discontinuity in the ecosystem response (Hillebrand et al. 2020). The increase of cumulative impacts can also create stressor interactions due to their spatial and temporal overlap (Côté et al. 2016). This situation can result in the emergence of complex relationships between stressors, known as synergistic (i.e., greater than the additive effect) or antagonistic interactions (i.e., less than the additive effect) (Piggott et al. 2015;Côté et al. 2016). Both threshold effects and stressor interactions create nonlinear biological responses, leading to "ecological surprises" in natural ecosystems.
Estuaries are among the ecosystems most affected by environmental stressors (Elliott and Quintino 2007). Already naturally stressed, with high variability in environmental factors, estuaries are highly dynamic ecosystems and are simultaneously exposed to high degrees of anthropogenic stress, including global changes (Elliott and Quintino 2007). Furthermore, human activities are often concentrated along the coasts, leading to higher cumulative exposure in these regions, with hotspots near coastal cities (Beauchesne et al. 2020). Among the current stressors in estuaries, particularly in shallow coastal areas, nutrient enrichment and subsequent eutrophication has gained much attention over the last decades (de Jonge et al. 2002;Alsterberg et al. 2012). Primarily due to the use of artificial fertilizer, nutrient enrichment may have direct or indirect impacts on organisms by generating anoxic/hypoxic zones or the occurrence of harmful algal blooms (de Jonge et al. 2002). These direct and indirect impacts can ultimately change the structure, the composition and the overall ecophysiological state of aquatic communities, such as benthic communities (de Almeida et al. 2007;Levin et al. 2009;Patil 2011), that can then alter fundamental ecosystem processes and services (Villnäs et al. 2012;Gammal et al. 2022). The influence of saltwater intrusion has also been a main concern for estuarine ecosystems (Velasco et al. 2019). The salinization of freshwater ecosystems, such as the fluvial portions of estuaries, can be caused by natural phenomena (e.g., weathering of the catchment, small amounts of salts dissolved in rainwater; Herczeg et al. 2001) and anthropogenic activities (e.g., irrigation for agriculture, sewage, and industrial effluents; Cañedo-Argüelles et al. 2013). In addition, sea level rise associated with climate change can also cause saltwater intrusion in freshwater ecosystems (Morrissey et al. 2014). Saltwater intrusion will have various effects on organisms/individuals (e.g., survival, fecundity, metabolic rates; Velasco et al. 2019), on population/community responses (e.g., predation, competition; Bäthe and Coring 2011;Mill an et al. 2011) and impact community structure (Tully et al. 2019) and ecosystem processes (e.g., breakdown of allochthonous organic matter; Fritz et al. 2010;Schäfer et al. 2012).
Many studies have characterized the interacting effect of multiple stressors and identified thresholds along increasing environmental gradients. While those two topics (i.e., stressors interactions and threshold effects) have been the scope of many experimental, modeling or observational work, very few studies have attempted to combine both and investigate stressor interactions along environmental gradients and thresholds in the context of single or multiple stressors (Carrier-Belleau et al. 2021b). Based on the environmental conditions of the St. Lawrence estuary (Quebec, Canada), this study aimed to determine the short-term effects of individual and combined stressors on multiple components of a simplified freshwater community in a laboratory setting. More specifically, our objective was to identify the effects of osmotic stress (i.e., saltwater input-pulse stressor) along a nutrient enrichment gradient (press stressor) that mimicked anthropogenically induced nutrient inputs on metabolic rates in zebra mussels, Dreissena polymorpha (Pallas 1771), on photosynthetic periphyton, and heterotrophic microbial activity. Finally, we aimed to identify breakpoints along the gradient in the context of single (i.e., stable salinity across different levels of nutrients) and multiple (i.e., variable salinity across different levels of nutrients) stressors. We hypothesized that saltwater input would interact with nutrient enrichment along the gradient, creating synergies or antagonistic interactions (Carrier-Belleau et al. 2021a). The effect would differ according to the response variable, knowing stressors interact differently depending on the level of biological organization (Galic et al. 2018). More specifically, we hypothesized nutrient enrichment would alter metabolic rates in D. polymorpha through direct (i.e., toxicity; de Almeida et al. 2007;Patil 2011) or indirect effects (i.e., lower dissolved oxygen concentrations; Burky 1983;Mcmahon 1996) and increase periphyton biomass (McDowell et al. 2020). On the contrary, as nutrient enrichment has been shown to both increase (Qualls 1984) and decrease microbial activity (Woodward et al. 2012;Kearns et al. 2016), we expected nutrient enrichment could have either positive or negative effects on decomposition rates in our experiment. Moreover, we hypothesized saltwater input would alter metabolic rates in response to stress (Kilgour et al. 1994), negatively affect phototrophic periphyton biomass through cellular disruption (Sudhir and Murthy 2004;Hu and Schmidhalter 2005) or shifts in assemblages (Mazzei et al. 2018) and decrease microbial activity as lower decomposition rates have been observed in natural and laboratory settings (Fang et al. 2010;Morrissey et al. 2014). Finally, we hypothesized that the breakpoint along the gradient might only be identifiable with the occurrence of osmotic stress (i.e., in the context of multiple stressors), as a growing body of evidence indicates that interactions between "pulse" stressors (i.e., temporary disturbances) and slow change "press" stressors (i.e., long term disturbances) could create breakpoints (Dudney and Suding 2020;Hillebrand et al. 2020). The overarching goal of this study is to demonstrate the importance of considering stressor interactions along gradients and threshold effects in the context of single and multiple stressors and their impact on multiple components of a freshwater ecosystem.

Choice of species, responses, and biological components
We investigated metabolic responses in D. polymorpha, changes in phototrophic periphyton biomass and microbial activity through an organic-matter decomposition assay. The choice of these taxa is based on the relative sedentary behavior which makes them exposed to an increasing amount of natural and anthropogenically induced stressors. Furthermore, we investigated changes in phototrophic periphyton biomass, as periphyton mats represent a unique microhabitat for many autotrophic and heterotrophic organisms and play an important role in food web dynamics, especially for benthic organisms (e.g., grazer, filter feeders; Hillebrand and Kahlert 2001). Periphyton biomass and productivity, particularly benthic algal assemblages, can be altered by environmental factors, such as salinity and major nutrient concentrations (i.e., nitrogen and phosphorus; Leland and Porter 2001). Finally, measuring microbial activity through an organic-matter decomposition assay (after Tiegs et al. 2019) allows the quantification of the capacity of the simplified community to process organic carbon when exposed to nutrient enrichment and saltwater input. Increased nutrient concentrations have become a central theme in freshwater ecosystems (Vitousek et al. 2009), and decay rates are known to vary along nutrient enrichment gradients, with slower rates at high and low nutrient concentrations (Woodward et al. 2012;Costello et al. 2022). Both phototrophic periphyton biomass and heterotrophic microbial activity were measured to assess the effect of combined stressors on ecosystem processes (Young et al. 2008).

Periphyton, microbial community, and specimen collection
Periphyton represent a major food source for macroinvertebrates, and periphyton biomass and productivity have been shown to vary along phosphorus (P) enrichment gradients in different freshwater ecosystems (McCormick et al. 2001). Initial periphyton and microbial communities were obtained using Hester-Dendy samplers (Hester and Dendy 1962). To mimic a natural periphyton and microbial community, we deployed 36 sets of 14 7.62-cm round masonite disks spaced 0.635 cm apart (Hester Dendy, NKY Environmental Supply, KT) in the Bassin Louise (46 31 0 N, 71 20 0 W) in Quebec City (Quebec, Canada) 4 months prior to the experiment (June 2020) to allow colonization. In mid-October, 1 month prior to the start of the experiment, plates were retrieved and taken to the laboratory at Université Laval (Quebec, Canada) where they were inspected to remove all macrofauna while retaining periphytic communities. Plates were placed in 37 L coolers which served as temporary aquaria to enable acclimation to laboratory conditions. Each aquarium was equipped with an air pump (Marina 50,11,110,Hagen) and a filter (Marina Slim Filter S10, A285) and was exposed to a 12 h photoperiod (07:00 h-19:00 h). Three plates were then randomly distributed in each experimental unit 1 week prior to the start of the experiment (see below).
Zebra mussels were collected using the same Hester-Dendy sampling plates used for periphyton colonization. Specimens were removed from the plates, and only individuals measuring less than 9 mm, for standardization, were kept for the experiment. A total of 3500 individuals were placed in the same coolers that served as aquaria during the acclimation period. Individuals of D. polymorpha were fed daily with Shellfish Diet 1800 (Reed mariculture, Campbell, CA). A week prior to the experiment, 1950 individuals were tagged on the right outer shell using bee-marking numbers (The Bee Works, Orillia, ON, Canada). All specimens were then transferred to the Laboratoire de Recherche en Sciences Aquatiques (LARSA; Université Laval) where we conducted the experiment and where we randomly distributed 15 mussels in each experimental unit (3.8-L plastic aquaria).

Experimental design and protocol
We performed a 2-month-long experiment between 16 November 2020 and 11 January 2021. To test for the effects of nutrient enrichment alone and coupled with salinity variation on responses at multiple levels of biological organization, we randomly assigned each experimental unit to one of the 26 experimental treatments following a full-factorial design (n = 5 replicates per treatment; Table 1). We used an in situ enrichment technique with controlled-release pellets in mesh bags (Acer NT NPK 17-7-10; Plant-Prod, Brampton, ON, Canada; Alsterberg et al. 2012;Carrier-Belleau et al. 2021a) to create a eutrophication gradient (13 levels) ( Table 1). The fertilizer pellets were placed in mesh bags attached from the surface of the aquarium, so the mesh bags all hung in the center of each aquarium. Nutrient loss was judged not to differ across treatments as fertilizer pellet percentage loss was roughly equivalent for all nutrient and salinity treatments (Table S1) and was similar to those used in other experiments simulating nutrient enrichments during the comparable time exposures (Baden et al. 2010;Carrier-Belleau et al. 2021a). The nutrient concentrations used mimicked a freshwater ecosystem passing from a mesotrophic to eutrophic state due to nutrient enrichment induced by agricultural activities (Smith et al. 1999). The freshwater reserve at LARSA is supplied by a municipal water intake and first passes through a sand filter (MAR-24-0, Culligan, QC, Canada) in order to eliminate all suspended particles and continues to a dechlorinator (MAR-C-42). Finally, a UV filter acts as a third filtration system to kill any pathogens that may have resisted the previous treatments. To identify how a second stressor can alter the presence of a tipping point, we exposed half of the aquariums to the same 13 levels of nutrient concentrations coupled with salinity variation. The stress induced by salinity variation was obtained by manually adding 400 mL of artificial saltwater to half of the aquariums once a week. The saltwater at LARSA is obtained by mixing mineral chemical and inorganic compounds to attain 30 PSU (Table S2). A mean salinity of 3 PSU was obtained in these experimental units and it took 4 h for the salinity to be re-established to initial conditions. This treatment was used to represent the shifting of an estuary salinity front, as often seen in estuaries (Liu et al. 2019). This salinity was also chosen in order to induce a second stress on D. polymorpha, while remaining in their range of tolerance (Wright et al. 1996). For control conditions (0 g of fertilizer pellets, stable salinity), we placed mesh bags filled with aquarium gravel to control for possible effects of mesh bags in aquaria and 400 mL of freshwater was added to the aquaria with stable salinity to control for possible effect of increased flow. Thus, we had a total of 26 experimental treatments, with five replicates per treatment for a total of 130 aquaria sustained by a constant water flow.
The 130 aquaria were randomly distributed in six water baths that maintained the temperature at 8 C. Temperature, salinity (during and in the absence of osmotic stress), pH (during and in the absence of osmotic stress), and dissolved oxygen (just before and during light exposure) were measured 3 days per week in all experimental units. Nutrient concentrations (nitrate and phosphate) were measured twice during the experiment (14 December and 8 January). To determine the nutrient concentrations, three water samples were collected in every aquarium 7.6 cm above the bottom of every aquarium. Water samples were filtered through GF/F filters (nominal porosity 0.7 μm) in an acid-rinsed falcon tube (15 mL) with an acidrinsed plastic syringe. Because of high concentrations, samples were then diluted (1 : 5) with MilliQ water and directly analyzed using a Bran and Luebbe Autoanalyzer 3 applying colorimetric methods adapted from Grasshoff et al. (1999). The mean values of experimental parameters are given in Table 1. Dead mussels were removed daily from the tanks to avoid increasing bacterial biomass in the aquaria. To keep the same abundance and density of mussels throughout the experiment, mussels from an extra set of aquaria for every treatment, each containing 100 untagged mussels, were used to replace dead individuals in the experimental aquaria with individuals exposed to the same conditions. These individuals were not, however, used for subsequent analyses.

Metabolic rate measurements of D. polymorpha
Oxygen consumption (as a proxy for metabolic rates) was measured using a closed-chamber respiratory method. The day before the end of the experiment, no food was added to the aquaria to starve mussels to prevent the effects of dynamic action on the metabolic rate measurement (Speakman and McQueenie 1996). Three individuals from each aquarium were then haphazardly collected, gently cleaned to remove any biofilm or periphyton, and carefully pooled in 6.3 mL borosilicate glass vials filled with water filtered using GF/F filters (nominal porosity 0.7 μm) to prevent bacterial and algae respiration and oxygen production (n = 5 replicates per treatment). Each vial was equipped with a small magnetic stir rod isolated from the organisms with a small mesh to ensure homogenous oxygen concentrations in the vials. The vials were placed in a water bath to maintain a constant temperature of 8 C and onto magnetic stirrer plates (Mix 15; 2Mag, AG, Munich, Germany) to activate stir rods. Oxygen concentration was measured using a noninvasive fiber-optic system (Fibox 4; PreSens, Regensburg, Germany), and reactive oxygen sensor spots were glued inside each vial. Linear calibration was previously achieved between the oxygen concentration in two solutions of bubbled overlying water (100%) and sodium ascorbate solution (0%), respectively. Measurements were recorded every 20 min and incubations were stopped after 30% oxygen consumption, around 70% O 2 saturation (typically after 5 h of incubation). The measurements were performed over 3 days, and five controls (vials containing only filtered water) were performed every day to control for potential microbial activity. Microbial oxygen uptake rates were very low (less than 2% over 5 h incubations) and were thus ignored in calculating oxygen uptake in D. polymorpha. Following the incubation period, shell metrics (length, height, width) and mass of each individual were taken. Individuals were then dissected, and the mass of each tissue was measured (blotted weight). The results are expressed as oxygen uptake in μmol total g of tissue À1 h À1 .

Phototrophic periphyton biomass
At the end of the experiment, the Hester-Dendy samples were retrieved from the aquaria and were scraped to collect the periphyton. For every plate, half of the surface was used for chlorophyll a (Chl a) and phaeopigment contents and to evaluate the ratio chla/phaeopigments (i.e., indicator of the physiological state of macroalgae). Periphyton samples (mean AE SE: 23.91 AE 0.5645 mg) were placed in 15 mL Falcon tubes and frozen at À80 C for subsequent analysis. The other half of the plate was preserved at À80 C to be used for other analyses not included in this study. Periphyton samples were placed in 10 mL of 90% acetone for 24 h at 4 C and the supernatant was analyzed fluorometrically (Turner 10 AU; San Jose, CA) before and after acidification (HCL 5%) (Riaux-Gobin and Klein 1993). After pigment determination, we used preweighed GF/F filters to collect the organic material (n = 15 per treatment) and determine the dry weight (mg g À1 ; Riaux-Gobin and Klein 1993; Link et al. 2011).

Heterotrophic microbial activity
To assess microbial activity, we measured organic-matter decomposition using a standardized cotton strip assay (Tiegs et al. 2013). Cotton strips are composed of greater than 95% cellulose, and the quantification of tensile strength loss of the cotton fabric reflects the microbial catabolism of cellulose, providing a proxy for organic-matter decomposition (Tiegs et al. 2013). The strips were prepared from bolts of Fredrix-brand unprimed 12-oz. heavy-weight cotton fabric and cut into 25-mm-wide (27-thread-wide), 8.0-cm-long cotton strips. Three strips were placed in each aquarium. To do so, a small hole was made with a thin nail 7 mm from the end of each cotton strip. A thin cable tie was then pushed through the hole and was attached to a wire crossing the width of the aquaria. The cable tie was thus in the air, whereas the whole strip was immersed in the water of the aquaria.
Strips were removed at the end of the experiment, and each strip was placed in 80% ethanol for 30 s before being dried at 40 C. Tensile strength of each strip was measured by following the protocol of Tiegs et al. (2013). Briefly, each strip was placed into the grips of a tensiometer (Mark 10 MG100) mounted to a test stand (Chatillon TCM 201) and pulled at a 20 mm min À1 rate. The peak tensile strength value was recorded for each strip.
The decay rate was calculated using the formula: where tensile strength treatment strip is the maximum tensile strength recorded for each strip incubated in our system (n = 15 per treatment) and tensile strength initial is the mean tensile strength of 10 strips that were not incubated (mean AE SD = 67.3 AE 0.9).

Statistical analysis
The effects of single and combined stressors were analyzed using linear mixed-effect models, where the terms "nutrient" (13 levels) and "salinity" (2 levels) were included as fixed factors and the terms "water bath" and "aquarium" were added as random factors with "aquarium" nested in "water bath." To test the effect of the random effects, we performed a variance component test with the log-likelihoods of the model with and without the random effects (varCompTest function of the varTestnlme package in R). Both random factors had no significant effect (p > 0.1). Response variables (i.e., metabolic rates in D. polymorpha, periphytic pigments and microbial activity) were first tested for the assumption of normal distribution with a Shapiro-Wilk test. Regression diagnostic plots were used to verify the homoscedasticity of residuals (equal variance of residuals). We set the possibility of making a type I error to α = 0.05 but considered p < 0.1 as marginal effects. We used pair-wise multiple comparisons when we obtained significant or marginally significant interactions among given sources of variations to understand which levels differed.
Breakpoints for all measured responses were determined using broken-line regression models with the SEGMENTED R package (Muggeo 2008; R Core team 2019). This type of analysis determined the relationship between observed responses and nutrient enrichment for both our stable and variable salinity treatments and identified changes in the regression relationship. The breakpoint along the gradient was then estimated by identifying a gap and a variation in the slope coefficient in the linear predictor (Vanacker et al. 2015). We also used the Davies test (Davies 1977) to determine if the slopes before and after the breakpoint were significantly different (p ≤ 0.05). Since breakpoint analysis is highly sensitive to data extremes, we identified and excluded outliers from the analysis using Dixon's test (Dixon 1950).

Interactions interpretation
To interpret the nature of the observed interaction between osmotic stress and different nutrient concentrations along our gradient, we used an additive (A + B) effects model (Piggott et al. 2015;Galic et al. 2018). This approach considers the deviation of the observed effect of combined stressors from a null model to determine interaction types (synergism and antagonism) and allows identification of the magnitude and direction (positive/negative synergism or antagonism) of the combined effect (Piggott et al. 2015;Jackson et al. 2016). To characterize the effect of nutrient enrichment at a given concentration and saltwater input, and to obtain our null additive model, we compared the observed response to our control conditions (no nutrient enrichment, stable salinity). The effect can be positive or negative: where CT is the observed response in control conditions, and N and SV are the observed effect (i.e., the difference between the response of a treatment and CT) of nutrient enrichment and salinity variation, respectively.

Metabolic rate
Oxygen uptake in D. polymorpha generally increased linearly (i.e., no breakpoints identified) between 0 and 110 g of fertilizer pellets in the context of stable salinity (Fig. 1), reaching a maximum at 110 g and slightly (but not significant) decreased past this point. In contrast, metabolic rate in the presence of osmotic stress reached two significant breakpoints along the gradient. A first breakpoint was identified at 30 g and a second at 100 g of fertilizer pellets. Thus, metabolic rate increased between 0 and 30 g, slightly decreased from 30 to 100 g and then decreased beyond this second threshold (Fig. 1).
Metabolic rates were also influenced by osmotic stress along the gradient, as shown by a significant interaction between nutrient enrichment and salinity variation (Table 2). Specifically, mean metabolic rate in D. polymorpha at 10 g of fertilizer pellets was greater for individuals exposed to osmotic stress, as shown by the significant difference between the stable and variable salinity treatments at that given nutrient concentration (pair-wise comparison, p = 0.0161). In fact, the observed combined effect of both stressors was greater than the effect of individual stressors and greater than the additive null model, indicating a positive synergistic interaction (i.e., less positive than the null model) between saltwater input and nutrient enrichment at 10 g of fertilizer pellets. Indeed, both nutrient enrichment and saltwater input reduced the metabolic rate compared to control conditions (Fig. 2a). However, in the presence of both stressors, oxygen uptake was seven times greater than the additive null model (0.0039 vs. 0.0005 μmol g À1 h À1 ) and twice the oxygen uptake in control conditions (0.0039 vs. 0.0019 μmol g À1 h À1 ). Mean concentrations for oxygen uptake are summarized in Table S3.

Microbial activity
The decay rate of cotton strips was less than 0.025 (tensile loss day À1 ) across all treatments (Fig. 3). There was no evidence of an influence of the nutrient enrichment gradient on microbial activity, but there was moderate evidence of an effect of osmotic stress (Table 2) with a lower microbial activity in the occurrence of variable salinity (Fig. 3). No breakpoints were identified along our gradient for both salinity conditions. Mean decay rates for all treatments are summarized in Table S3.

Phototrophic periphyton biomass
For all treatments, Chl a concentration was, on average, less than 9 mg per g of dry mass (Fig. 4a). It was negatively Fig. 1. Individual and combined effects of saltwater input along a nutrient enrichment gradient (initial g of fertilizer pellets) on oxygen uptake in Dreissena polymorpha. Points represent raw data for oxygen uptake observed in individual tanks. Full lines represent the broken-line regression models for the stable and variable salinity and shaded zones represent 95% confidence intervals. Dashed vertical lines represent significant (p < 0.05) thresholds along the gradient. No significant threshold was found for the stable salinity. Table 2. Results of linear mixed-effect models investigating the effect of nutrient enrichment (N + ) and salinity variation (S) on periphytic pigments, oxygen uptake in Dreissena polymorpha, and microbial activity. Details are provided for the F ratio and probability level (p value). affected by osmotic stress, as shown by a marginal effect of saltwater input (Table 2). On average, Chl a concentration was less for a given nutrient for the variable salinity (7.75 mg g À1 of dry mass) compared to the stable salinity treatment (9.67 mg g À1 of dry mass). This trend can be observed along the gradient for nutrient concentrations ranging from 30 to 80 g and 90 to 120 g of fertilizing pellets. Phaeopigments varied significantly as a function of the nutrient Â salinity variation interaction ( Fig. 4b; Table 2). Pair-wise comparison showed a difference between the two salinity treatments at 70 and 120 g (p = 0.0600 and 0.0468, respectively). At the two nutrient concentrations, nutrients and saltwater input both decreased phaeopigments concentrations (Fig. 2b,c). The combined effect of saltwater input and nutrient enrichment was greater than an additive null model, but lower than the effect of single stressors, resulting in a negative antagonistic interaction (i.e., less negative than predicted additively).

Source of variation
The ratio Chl a/phaeopigments was not affected by single or combined stressors along the nutrient enrichment gradient ( Fig. 4c; Table 2). Finally, no breakpoints were found for Chl a, phaeopigments, and the ratio Chl a/phaeopigments for both salinity treatments (Fig. 4c). Mean pigment concentrations for all treatments are summarized in Table S4.

Discussion
Identifying the effects of individual and combined stressors on multiple biological components and characterizing potential interactions along environmental gradients is a major challenge in natural ecosystems studies. Here, we aimed to identify stressor interactions along an environmental gradient and to characterize how multiple stressors can alter or create thresholds. We showed that thresholds and interactive effects are not ubiquitous and depend on the biological response and component investigated. However, we also showed that multiple stressors might generate threshold effects that complicate management actions. Thus, our results highlight the importance of considering both concepts (i.e., interactions among multiple stressors and threshold effects) in the context of multiple stressors to assess how interactions may either accentuate or create thresholds in natural systems. Below, we focus our discussion on the effect of individual or combined stressors depending on the biological compartment, the presence of thresholds along the nutrient enrichment gradient, and recommendations for future work. For the combined effect (N + Â S), a red bar corresponds with a positive synergy (S + ; more positive than predicted additively) and a blue bar represents negative antagonistic interactions (A À ; less negative than predicted additively). Error bars correspond to standard errors.

Effect of individual stressors
While stressors may interact, they may also have pronounced individual effects resulting in the dominance of a single stressor. For instance, as expected regarding microbial activity, we saw moderate evidence of a dominant negative effect of saltwater input, with reduced organic-matter decay rates in the presence of osmotic stress. As far as microbial activity is concerned, changes in salinity may impact extracellular enzyme activity by influencing molecular stability and protein confirmation states (Morrissey et al. 2014). Laboratory studies have shown a detrimental effect of salt input on enzyme activities (Das et al. 1997;Fang et al. 2010). The extracellular enzyme-mediated hydrolysis of substrates into monomers and oligomers represents a major aspect of microbial activity (Shi et al. 2011) and regulates decay rates in various habitats (Freeman et al. 2004;Allison and Vitousek 2005). This lower microbial activity, resulting in lower decay rates and organic-matter processing, can have major implications for estuarine systems by altering local conditions (e.g., nutrient cycling [P, N], fluxes of particulate and dissolved organic matter and primary productivity; Simon et al. 2002, Statham 2012. As hypothesized, a negative effect of salinity variation on Chl a concentration was also identified. This can be explained by the fact that saltwater intrusion can disrupt cellular ion homeostasis and result in decreased chlorophyll content, photosynthetic electron transport rate, and the uptake of essential nutrients (Sudhir and Murthy 2004;Hu and Schmidhalter 2005).
Changes in species composition, which were not measured here, may also explain the decrease in Chl a content and a significant shift in periphytic diatom assemblages have been reported due to saltwater intrusion (Mazzei et al. 2018). This shift in composition and relative abundance of diatoms or other mat-engineering species, such as cyanobacteria, can have significant implications by altering the periphyton mat type and, thus, the provided ecological service (Hagerthey et al. 2011).
In this case, dominant effects may still generate ecological surprises, as we expected to see a dominance of nutrient input on microbial activity (Qualls 1984;Kearns et al. 2016) and Chl a (Hao et al. 2020;Hope et al. 2020) but no effect was detected. First, in other studies, nutrient enrichment is associated with both positive and negative effects on microbial activity. For instance, higher decomposition rates have been associated with nutrient enrichment (Qualls 1984) and could be attributed to the microbial immobilization of exogenous inorganic nitrogen to build cellular proteins and hexosamines (Kaushik and Hynes 1971). On the other hand, Kearns et al. (2016) have shown that nitrogen additions can induce diversity loss in the active portion of the community through the increase in the proportion of dormant microbial taxa (Kearns et al. 2016). In our study, the lack of nutrients' effect on microbial activity could be explained by either/or (1) the fact that the range of tested nutrient gradient is above or below an optimal window of nutrients concentrations for microbial activity, (2) the algae present in the experimental unit may have competed with microbial heterotrophs for available Fig. 3. Individual and combined effects of saltwater input along a nutrient enrichment gradient (initial g of fertilizer pellets) on cotton strip decay rates, used as a proxy for microbial activity. Points represent raw data for decay rates observed in individual tanks. Full lines represent smoothing lines for the stable and variable salinity and shaded zones represent 95% confidence intervals. No thresholds were found for either salinity treatment.
nutrients in the mesocosms, and (3) factors other than nutrients modulated microbial activity (Webster and Benfield 1986;Woodward et al. 2012). However, at this time we do not have empirical evidence to support or exclude any of these potential mechanisms, and further studies will have to be designed to shed light on these. Second, we initially expected nutrient Fig. 4. Individual and combined effects of saltwater input along a nutrient enrichment gradient (initial g of fertilizer pellets) on (a) chlorophyll a concentration (b) phaeopigments concentrations and (c) the ratio chlorophyll a/phaeopigments. Points represent raw data for pigments (chlorophyll a or phaeopigments) observed in individual tanks. Full lines represent smoothing lines for the stable and variable salinity, and shaded zones represent 95% confidence intervals. No thresholds were found for either salinity treatment. enrichment to stimulate periphyton biomass, as it has been shown in laboratory and natural settings (Pinckney et al. 1995). However, the lack of nutrient effect on periphyton Chl a has been noticed in coastal aquatic systems, such as streams, where other environmental drivers (e.g., light) are the primary limiting factor (Kim and Richardson 1999). Among other things, this could also be explained by the presence of micrograzers in our mesocosms (Vermaat 1994).

Interactions among stressors
Stressors may interact and create synergistic or antagonistic relationships, depending on the biological compartment or response of interest. In our study, this phenomenon was first illustrated by a synergistic interaction between salinity variation and nutrient enrichment at 10 g of fertilizing pellets for the metabolic rate in D. polymorpha (Fig. 2a), where the oxygen uptake in the presence of both stressors was twice the uptake than in control conditions (no nutrients, stable salinity). At this nutrient concentration, both nutrient enrichment and salinity variation individually decreased oxygen consumption rates in comparison to control conditions. This was to be expected, as both stressors are known to individually alter metabolic rates in invertebrates. For instance, shell valve closure in bivalves has been proposed as a protective response to cope with osmotic stress or pollutants, allowing the organisms to survive under varying environmental conditions (Fournier et al. 2004;Liao et al. 2009). However, besides compromising food intake, valve closure also restricts water circulation and gas exchange and will progressively lead to the absence of oxygen availability in the mantle cavity, altering metabolism (Byrne et al. 1990). Ultimately, a lack of water circulation will also prevent the ejection of metabolic CO 2 and the maintenance of low osmolality of body fluids through ions transport and acid-base balance, which will have an impact on metabolism (Byrne and Dietz 1997;Ortmann and Grieshaber 2003). To deal with environmental stress, mussels have also developed a series of adaptations allowing them to survive under stressful conditions, including changes to respiratory and overall metabolic rate, activation of alternative pathways for energy production, induction of biochemical defense and repair mechanisms (de Almeida et al. 2007).
Thus, individuals could be able to cope with individual stressors, but increase their metabolic rate significantly, indicating a potential response to stress, when exposed to a combination of stressors. For example, at this nutrient concentration (10 g) in the presence of osmotic stress, individuals of D. polymorpha could have increased their oxygen uptake, as bivalves are also known to increase their metabolic activity to survive extended periods of stress from pollutants (de Almeida et al. 2007;Patil 2011). This can be explained by the fact that pollutant stress will cause high energy expenditure for active elimination or detoxification (de Zwaan et al. 1995). We cannot exclude, however, that an increase in the metabolic rate at 10 g of fertilizer pellets increased organismic aerobic scope and/or its net ATP production (Bozinovic 2011). Thus, the synergy observed at this nutrient concentration may imply a beneficial impact on individuals of D. polymorpha. The absence of synergy at other levels of nutrient enrichment make it even more challenging to (1) find a plausible mechanistic explanation for the synergy observed at 10 g of fertilizer pellets and (2) to conclude if an increased metabolic rate in this situation translates into a cost for the organisms or on the contrary, in an advantageous effects. Further studies along a broader nutrient enrichment gradient could allow us to identify other interactions (synergies or antagonisms) along the gradient which would facilitate the mechanistic interpretation. On the other hand, investigating other responses, such as energy content, filtration rates or reproduction for example could allow us to interpret metabolic rates results.
Antagonistic interactions between salinity variation and nutrient enrichment were identified at 70 and 120 g of fertilizer pellets for phaeopigment concentrations. In fact, at both nutrient concentrations, both stressors, alone and in combination, reduced phaeopigment concentrations, and, thus, Chl a degradation. At these nutrient concentrations, nutrient enrichment could have had beneficial effects by decreasing the overall degradation of chlorophyll pigments (Bricker et al. 2008). On the other hand, salinity could have decreased meiofauna (e.g., rotifers, nematodes, oligochaetes-Majdi et al. 2020) grazing, which were most certainly present in our communities, as salinity is a major environmental factor influencing the presence of benthic organisms (Attrill 2002), especially meiofauna (Zeppilli et al. 2018). When combined, stressors could act on these two mechanisms or another mechanism that we could not identify with our results.

Presence of thresholds and implications for management actions
The occurrence of breakpoints along our environmental gradient was only observed for metabolic rates in D. polymorpha under multiple stressor conditions, suggesting, in this case, acting upon saltwater intrusion might result in the highest ecological benefit. As explained above, different mechanisms could explain the presence of thresholds at 30 and 100 g of fertilizer pellets under osmotic stress, and their interpretation becomes even more difficult, as an increase or decrease in metabolic rates can both represent a stress response (Bozinovic 2011). However, we wish to emphasize that our results show that the presence of multiple stressors (i.e., nutrient enrichment and saltwater intrusion) changes the relationship between a biological response and an environmental gradient from being linear to experiencing threshold effects. Biological responses tend to be more linear and question the application of thresholds in natural systems (Hillebrand et al. 2020). Others propose that thresholds should not be expected along singular gradients but could be expected when organisms or communities are exposed to the combination of pulse "triggers" and slowly changing stressors or environmental gradients (Dudney and Suding 2020). Threshold effects have long been used in ecosystem management, conservation, and restoration (Vaquer-Sunyer and Duarte 2008; Suding and Hobbs 2009), as nonlinear effects are harder to predict and understand. However, our results show that we now need to integrate threshold effects in the context of multiple stressors to guide these actions. For example, salinity intrusion, here treated as a pulse stressor, may have huge impacts on various biological responses and components in estuarine systems, as observed in this study. Salinity intrusion has a stronger negative individual effect on organismal performance traits than other stressors, such as temperature or toxicants (Velasco et al. 2019), highlighting the importance of increasing management efforts for this single stressor. Knowing the underlying causes of the temporal and spatial dynamics of this stressor in ecosystems may help establish precautionary principles for environmental policy and anticipate major transitions in natural systems (Vaquer-Sunyer and Duarte 2008; Suding and Hobbs 2009).

Concluding remarks and recommendation for future work
Studying stressor interactions along environmental gradients is a much-needed step toward integrating the concepts of stressor interactions and threshold effects. Experimental approaches, such as the work done here, may help identify (1) the dominance of single stressors, (2) interactions among multiple stressors, and (3) threshold effects in the context of single and multiple stressors. These components may, afterwards, be integrated into decision tools in local and global management.
From this experiment, some considerations for future research arise. First, work should concentrate on multiple stressors, as ecosystems are now mostly affected by cumulative impacts (Halpern et al. 2019). This can be done by focusing on stressors that are known to have dominant effects, interact with others, or create threshold effects. For example, our work highlights the dominance of saltwater intrusion on multiple ecological components and biological responses, suggesting that acting upon this dominant stressor at local and regional scales would result in the largest benefit for the ecosystem and eliminate potential threshold effects. In order to guide management actions, multiple-stressor research can be done by performing studies (e.g., laboratory or in situ experiments, observational studies, modeling) that include both local and global stressors, which require management at vastly different scales but are necessary to consider in tandem. Second, research has focused on identifying thresholds along singular gradients, but we recommend concentrating efforts on combining the ecological concepts of thresholds with the concepts of stressor interactions. This will allow identification of different interactions for stressor pairs and the switch between antagonistic to synergistic, or vice versa (Jackson et al. 2016). This will strengthen the impact of experimental studies because it better reflects conditions in natural ecosystems. Finally, stressors can affect different pathways of actions and differentially affect the responses measured at the individual, population, or ecosystem level (Galic et al. 2018) or when considering different species, trophic groups, or habitats (Gissi et al. 2021). The effect of combined stressors and the identification of thresholds should thus be considered for different biological compartments (i.e., from cellular responses to ecosystem processes). To do this and to maximize research efforts to guide management actions, studies should focus on specific species (e.g., commercial species, ecosystem engineers) but also focus on processes that are known to affect those species of interest (e.g., primary production, organic-matter remineralization). Combining multiple responses and components will allow us to adopt a holistic vision and approach to tackle the effect of multiple stressors in aquatic ecosystems. This will provide useful information for managers to identify the component for which action will need to be prioritized.