Stressors Increase the Impacts of Coastal Macrofauna Biodiversity Loss on Ecosystem Multifunctionality

There is substantial evidence that biodiversity underpins ecosystem functioning, but it is unclear how these relationships change with multiple stressors in complex real-world settings. Coastal zones are affected by numerous stressors (for example, sediment input and nutrient runoff from land) and the cumulative effects of these stressors may result in pronounced and unexpected changes in the functioning of ecosystems. To investigate the cumulative effects of turbidity and elevated nutrients on coastal biodiversity-ecosystem functioning relationships, we performed a large-scale field experiment manipulating in situ sediment porewater ammonium concentrations and measured multiple ecosystem functions related to carbon fixation and mineralisation in 15 estuaries with varying levels of turbidity. The results indicated that the benthic macrofauna diversity (species richness, abundance, and functional richness) declined with increased porewater ammonium concentrations and there were clear thresholds in light levels at the seafloor in relation to the biodiversity-ecosystem function relationships. Multifunctionality indices (an integrated index of all measured functions) in moderately turbid and turbid estuaries (daily mean seafloor PAR < 420 µmol m−2 s−1) decreased with the loss of macrofauna biodiversity. Functioning in low-turbidity estuaries (daily mean PAR > 420 µmol m−2 s−1) however remained relatively constant, suggesting that they were more resilient against the nutrient-induced biodiversity loss. Our results demonstrate that ecosystems already stressed by stressors that alter functional performance (turbidity) may be more prone to loss of overall functioning if biodiversity is reduced by another stressor (nutrient enrichment), highlighting the potential snowballing effects of cumulative change.


ABSTRACT
There is substantial evidence that biodiversity underpins ecosystem functioning, but it is unclear how these relationships change with multiple stressors in complex real-world settings. Coastal zones are affected by numerous stressors (for example, sediment input and nutrient runoff from land) and the cumulative effects of these stressors may result in pronounced and unexpected changes in the functioning of ecosystems. To investigate the cumulative effects of turbidity and elevated nutrients on coastal biodiversity-ecosystem functioning relationships, we performed a large-scale field experiment manipulating in situ sediment porewater ammonium concentrations and measured multiple ecosystem functions related to carbon fixation and mineralisation in 15 estuaries with varying levels of turbidity. The results indicated that the benthic macrofauna diversity (species richness, abundance, and functional richness) declined with increased porewater ammonium concentrations and there were clear thresholds in light levels at the seafloor in relation to the biodiversity-ecosystem function relationships. Multifunctionality indices (an integrated index of all measured functions) in moderately turbid and turbid estuaries (daily mean seafloor PAR < 420 lmol m -2 s -1 ) decreased with the loss of macrofauna biodiversity. Functioning in low-turbidity estuaries (daily mean PAR > 420 lmol m -2 s -1 ) however remained relatively constant, suggesting that they were more resilient against the nutrient-induced biodiversity loss. Our results demonstrate that ecosystems already stressed by stressors that alter functional performance (turbidity) may be more prone to loss of overall functioning if biodiversity is reduced by another stressor (nutrient enrichment), highlighting the potential snowballing effects of cumulative change.
In low-turbidity estuaries, benthic functioning was resilient to biodiversity loss Impacts of biodiversity loss increase in estuaries experiencing multiple stressors

INTRODUCTION
Alarming rates of global biodiversity loss have increased the need to assess the relationship between biodiversity and ecosystem function and service delivery (Cardinale and others 2012;Balvanera and others 2014;IPBS 2019). Knowledge regarding biodiversity-ecosystem function relationships has significantly increased during recent decades through various experiments that have included a range of taxa, trophic levels, habitats, and ecosystem functions (for example, Balvanera and others 2006;Gamfeldt and others 2015;Lefcheck and others 2015; van der Plas 2019; Lam-Gordillo and others 2020). There is thus a substantial amount of evidence that biodiversity is important for ecosystem functioning, but much of the research has been limited to small-scale experiments with low species diversity, single functions, and constrained environmental conditions (that is, in a controlled laboratory setting or one place and time; Snelgrove and others 2014). Such constrained experiments are highly context dependent and lacking generality with uncertainties in how these relationships change along spatial and temporal environmental gradients and particularly how they respond to multiple stressors (Thrush and Lohrer 2012;Brooks and Crowe 2019;O'Brien and others 2019). In times of rapid environmental change and a biodiversity crisis, there is a need to further resolve the types of biodiversity-function relationships that are generated by a multitude of complex connections between diverse species and functions that exist in natural ecosystems.
Often the relationship between cumulative stressors and ecological responses is not well resolved, but at the heart of responses are the interactions between taxa, processes, and the abiotic environment. Further complexity is added by the fact that ecosystems simultaneously provide many functions (that is, multifunctioning; Byrnes and others 2014; Manning and others 2018). Thus, studying biodiversity-multifunctioning relationships in ecosystems with diverse functions, species and connections informs a more general understanding of how ecosystems respond and can be resilient to stress (Siwicka and others 2020;Siwicka and others 2021). This type of ecological knowledge is needed to inform the potential for cumu-lative effects of multiple stressors, information that is critical to environmental management (Hewitt and Thrush 2019). Despite the important role of biodiversity and ecosystem function in understanding resilience to stress, there are few examples that link biodiversity to multiple functions and explore how these relationships change with stress in real-world settings (but see for example, Brooks and Crowe 2019). Here, we analysed data from a large-scale coastal marine field experiment that measured intertidal soft-sediment species diversity and multiple functions across two stressor gradients, turbidity and nutrient enrichment. Our study aimed to explore how biodiversity-ecosystem function relationships change over gradients of stress.
Coastal ecosystems are heavily impacted by many stressors originating from land-based activities, such as agriculture, horticulture, urban development, and forestry (Halpern and others 2015). Nutrient and sediment loading into coastal ecosystems have emerged as major drivers of coastal ecosystem degradation (Nixon 1995;Levin and others 2001;Valiela and Bowen 2002). These two stressors have multiple ecological effects that generate complex responses when they co-occur. Sediment runoff results in increased turbidity in the water column, and the turbidity hinders light from reaching the seabed which alters benthic primary production and thus have cascading effects on ecosystem carbon budgets and food webs (Miller and others 1996;Duarte and others 2005; Christianen and others 2017). Increasing nutrient additions in low concentrations can have positive effects on primary production, particularly in oligotrophic nutrient limited waters, but at higher concentrations, nutrients can cause eutrophication, changes to organic matter quality and quantity, hypoxia and ammonia toxicity which will change the functioning of the ecosystems as well as alter the macrofaunal communities (Kohn and others 1994;Douglas and others 2017;Thrush and others 2017;Breitburg and others 2018). The multiple effects of these two stressors on both the biota and the biogeochemical cycles result in changes to the linkages between species, environmental properties, and processes (biogeochemical and primary production).
The benthic macrofauna play a central role in regulating how the ecosystem functions as many of them alter sediment biogeochemistry and the exchange of solutes and particles across the sediment water interface (Kristensen and others 2012;Volkenborn and others 2012;Vanni and McIntyre 2016). The specific functional role that each species has on the ecosystem depends on its functional traits: location in the sediment (surface or deep), its movement (upwards, downwards, sideways), and feeding behaviour (for example, suspension or deposit feeding). Bioturbation, by macrofauna that mix or irrigate the sediment with their movements, enhances microbially driven remineralisation processes through oxygenation of the sediment as well as transport of particles and solutes within the sediment and across the sediment-water interface (Glud 2008;Kristensen and others 2012). Their activities are thus tightly linked to enhanced benthic primary production in oligotrophic systems through the release of porewater nutrients from deeper sediments (Lohrer and others 2005;Jones and others 2011;Lohrer and others 2015).
The multitude of different macrofauna species, that make up benthic community diversity, has different effects on their surroundings and functioning, which means that context dependencies in stressor impacts on function are likely to vary depending on the resident community and which species are lost (Gladstone-Gallagher and others 2019; Hewitt and others 2019). A diverse community is more likely to have functional redundancy (that is, species that can substitute for each other), and the communities' response diversity is also important for promoting functional resilience (that is, diversity of responses among species that contribute to the same function) (for example , Walker 1992;Yachi and Loreau 1999;Mori and others 2013). Concepts like response diversity and functional redundancy that occur in diverse communities complicate biodiversity-ecosystem function relationships, because in the real-world species and corresponding functional losses are nonlinear (Hagan and others 2021).
Despite an impetus to develop biodiversityecosystem functioning studies from smaller to larger scales within real-world settings, empirical knowledge of multiple stressor effects on the relationships remains scarce especially in marine ecosystems (Gamfeldt and others 2015; Lefcheck and others 2015; Brooks and Crowe 2019; O'Brien and others 2019). Moreover, there is also a need to move away from measuring impacts of biodiversity loss on single functions and use methods that incorporate the complexity of ecosystems that are underpinned by many interacting functions (Byrnes and others 2014; Dooley and others 2015; Manning and others 2018; Siwicka and others 2020; Siwicka and others 2021). One way of approaching this challenge is to utilise multifunctionality indices to link univariate measures of functioning to biodiversity measures (Hooper and Vitousek 1998; Byrnes and others 2014) within different environmental conditions. Manipulative experiments embedded in natural heterogeneous habitats can be valuable to advance our understanding of ecosystem function and multiple stressor effects (Bracken and others 2008;Snelgrove and others 2014;Thrush and others 2021).
To advance understanding of how biodiversity ecosystem function relationships are altered by stressor interactions, we undertook a large-scale field experiment. Porewater nitrogen concentrations were manipulated at 23 sites in 15 estuaries across New Zealand with different availability of light at the seafloor to quantify effects on biodiversity and ecosystem multifunctionality. The sites were carefully chosen to encompass a wide gradient of turbidity, but also to only include sites where populations of large shellfish occurred, indicating that sites were in similar habitats and not severely degraded. With this design and site selection of systems differently affected by turbidity, we aimed to specifically examine the ecosystem responses when adding a second stressor, nutrient enrichment, and subsequently to generalise the patterns of cumulative effects of multiple stressors on an intertidal ecosystem before the system has totally collapsed. We expected that the form of the biodiversity-ecosystem function relationships would change with the stressor gradients, but due to the heterogeneous environmental settings and potential complex interactions, the forms of those relationships could not be hypothesised a priori.

Study Design
Sediment runoff from land is a large threat for coastal ecosystems in New Zealand (Morrison and others 2009;MacDiarmid and others 2012) and the turbidity in the water column is mostly caused by sediment input (a function of land use, soil type, climate and topography) and further resuspension of previously received sediment (for example, Vant and others 1990). Across the length of New Zealand, 23 sites in 15 estuaries encompassing a variation in water column turbidity (range of daily average seafloor light intensity 134-712 lmol photons m -2 s -1 ) were identified ( Figure 1). The study estuaries were all relatively oligotrophic with low water column nutrient concentrations (at the time of sampling NH 4 + 0.04-2.8 lmol l -1 , NO x and P mostly below detection limits). The sites were located in sandy mid-intertidal sediments with natural densities of the bivalves Macomona liliana Stressors Increase the Impacts of Coastal Macrofauna… and Austrovenus stutchburyi. The sites were specifically chosen so that these functionally important species were present, indicating sites are in similar habitats and not heavily degraded (Clark and others 2020).
At each site, porewater nitrogen concentration was manipulated in nine 9-m 2 plots. Three plots were allocated to each treatment: medium nutrient addition 150 g N/m 2 , high nutrient addition 600 g N/m 2 , and procedural controls (no nutrient addition). In the nutrient addition plots, nitrogen (urea) was added uniformly into the sediment through injection of a slow-release fertiliser (40-0-0 N:P:K, Nutricote Chisso-Asahi Fertilizer Co., Ltd, Tokyo, Japan). The injection of the fertiliser at the depth of 15 cm elevates the porewater ammonium concentrations for extended periods, thus simulat-ing the effects of porewater enrichment due to long-term eutrophication but not necessarily the processes of increased organic matter decomposition (Douglas and others 2016; Douglas and others 2018). The nutrient enrichment treatments modified benthic macrofauna communities without driving them to extinction (see Results and Douglas and others 2018). All sites and experimental plots were established in March-April 2017 and left undisturbed for seven months before sampling sediment properties, macrofaunal diversity metrics and ecosystem functions in October-November 2017, which is springtime in New Zealand. The experimental plots were left for seven months to provide time for recovery from the initial set up, for the added nitrogen to diffuse through the sediment and for the ecosystem to respond. Further details of this study's design and sampling can be found in Thrush and others (2021).

Environmental Variables
Photosynthetically active radiation (PAR) was logged at each site for the duration of the experiment to quantify seafloor light intensity. The attenuation of light through the water column to the seafloor is affected by water column turbidity, and seabed PAR levels affect rates of photosynthesis by benthic primary producers. OdysseyÒ Light loggers measured PAR every 10 min at 10 cm above the sediment surface for 8 months. After the experiment, PAR data during daylight periods of tidal immersion (± 2 h of high tide) were summarised into daily averages, and then, these daily averages averaged for each site (hereafter called mPAR); data presented in Mangan and others 2020). The nights were excluded through removal of times when PAR was 0 lmol photons m -2 s -1 .
Sediment characteristics were sampled in each plot by pooling five sediment cores (2.6 cm diameter, 2 cm deep) for analysis of sediment chlorophyll a, grain size, and organic content. An additional four sediment cores (2.6 cm diameter) were taken from each plot for analysis of porewater ammonium concentration at two depths (0-2 cm and 5-7 cm). These depths were selected to quantify the nutrient concentrations above and below the redox potential discontinuity layer, and to ensure that the added nutrients had diffused upwards through the sediment. Sediment was frozen and later analysed for sediment chlorophyll a, grain size, organic content and porewater ammonium concentrations using standard techniques for sediment analysis in intertidal soft-sediments (see Thrush and others 2021 for details). We focused on measuring the porewater concentration of ammonium to quantify the treatment effect (that is, urea hydrolyses to ammonium), nitrate, nitrite, and phosphorous concentrations were however also measured in the porewater but not reported here due to being very low and mostly under the detection limits.

Biodiversity Measures
From each plot, two macrofauna cores (13 cm diameter, 15 cm deep, total area 0.026 m 2 ) were pooled, sieved on a 500-lm mesh and preserved in 70% isopropyl alcohol. The fauna were stained with Rose Bengal, sorted and identified. The deeper depth of the macrofauna cores, compared to the sediment and porewater samples, was used to ensure that all individuals were captured, the majority of the benthic macrofauna are however living in and affecting the sediment surface layers. To link biodiversity to the measured ecosystem functions, we calculated total number of species (S) and total abundance (N, number of individuals) per two cores (0.026 m 2 ) per experimental plot, as well as a measure of functional richness (F Ric ; Villé ger and others 2008; Laliberté and Legendre 2010). F Ric calculations were based on known biological traits of the taxa, which are expected to affect remineralisation processes in the sediments by pumping pore water, moving sediment particles and organic matter, and changing sediment topography (Villnä s and others 2012; Volkenborn and others 2012; Woodin and others 2016). These traits included categories of living position, sediment topographic features created, direction of sediment particle movement, motility, feeding behaviour, body size, and body shape (see Appendix 1 for traits and modalities used). As all the traits were numeric (probabilities with values ranging from 0 to 1), they were standardised to mean of 0 and unit variance and functional diversity was run based on Euclidean distances. Dimensionality reduction occurred with only the first 10 axes retained. F Ric output values were not standardised by the 'global' F Ric to run between 0 and 1. The F Ric calculations were performed using the dbFD package in R (Laliberté and others 2014).

Ecosystem Functions
We measured the rates of several processes that we use as proxies for ecosystem functions in each experimental plot (9 plots per site). Oxygen fluxes across the sediment-water interface were measured in situ using benthic chamber incubations. All incubations were performed on sunny days with mid-day high tides, to accommodate sampling at all sites this occurred over a five week period (26th October-27th November 2017) 7 months after the experiment was set up. In each experimental plot, two square metal frames (50 9 50 9 15 cm height) were pressed down 5 cm into the sediment at low tide. During the incoming tide, at water depths of approximately 50 cm, the chambers were sealed with transparent acrylic domes fixed to the metal frames, enclosing about 40 l of ambient seawater overlying the sediment. One chamber was covered with an opaque shade cloth in order to prevent light from entering the chamber (dark treatments) and the other was left uncovered to allow light to reach the seabed (light treatments). The chambers were incubated for approximately 4 h over a midday high tide. Water samples (1 9 60 ml) were collected through sampling ports at the start and the end of each incubation so that we could calculate a change in oxygen concentration in the water over 4 h (that is, a flux into or out of the sediment through time). The oxygen concentration in each water sample was measured with an optical probe (ProODO YSI, Yellow Springs, Ohio, USA), and the fluxes were calculated as (C end -C initial 9 V)/ A 9 t, where C is the oxygen concentration (lmol l -1 ), V is the volume of seawater inside the chamber (l), A is the area of sediment enclosed by the chamber (m 2 ), and t is the incubation time between initial and final samplings (h). The oxygen fluxes measured in the darkened chamber represent benthic respiration without photosynthesis, whereas the resulting fluxes from the light chambers represent net primary production (NPP) through photosynthesis. For further analyses with ecosystem functions, the dark chamber oxygen flux (that is, the respiration) was inverted to positive values and hereafter called sediment oxygen consumption (SOC). Gross primary production (GPP) was estimated as the sum of NPP and SOC in spatially paired chambers (paired by plot).
To measure rates of organic matter degradation at different sediment depths in the experimental plots, we used rapid organic matter assays (ROMA technique; O'Meara and others 2018). These assays involve the preparation and deployment of acrylic plates with 0.9 ml wells arranged in three vertical rows. The plates are inserted into the sediment with the wells corresponding to sediment depths of 1, 3, 7, 10, 15 cm into the sediment. Each well was filled with a 0.026 g C ml -1 mix of food grade agar, microcrystalline cellulose and powdered bran flakes. Carbon degradation rate was measured by the change in agar-mix volume in each well over time. One plate was deployed in each plot for the 5-8 d prior to the benthic chamber incubations. As a measure of carbon degradation at each site, we used the average carbon degradation rate (g C m -2 day -1 ) across 0-3 cm surface sediment depth (hereafter called CD). The values from the top 3 cm were used because within the sediment surface layers the highest rates of mineralisation processes are observed and most of the labile carbon is still considered to be in the system, as well as the most active layer with large quantities of benthic macrofauna (Chen and others 2017; Morys and others 2017).

Calculations of Multifunctionality
In this study, we quantified four processes that represent proxies for important ecosystem functions relating to carbon fixation and mineralisation. The four quantified processes were gross primary production (GPP), sediment oxygen consumption (SOC), carbon degradation (CD) and sediment chlorophyll a/phaeophytin ratio (Chl a/phaeo). GPP and SOC represent the uptake of CO 2 for primary production and total benthic respiration rates (an integrative measure of many processes (Glud, 2008)) across short time scales (that is, measured during a 4-h incubation). CD encompasses carbon consumption by both infauna and microbially driven mineralisation processes, within the surface sediment measured over about a 7-d period. Chl a/ phaeo in the surface sediment is used as a proxy for processes linked to autotrophic organic matter turnover (that is, the 'freshness' of the autotrophic organic matter, which is affected by both grazing and senescence; Bianchi and others 1988;Villanueva and Hastings 2000).
We used two methods to combine these four proxies and create indices that describe multifunctionality: (i) an average function approach (MF index ) and (ii) a thresholds approach (Gamfeldt and others 2008; Maestre and others 2012a; Maestre and others 2012b; Byrnes and others 2014). For the averaging approach, all functions were standardised by dividing by the maximum observed value, and the index (MF index ) was calculated by taking the average across all functions within one sample. The maximum value was based on the average of the five highest observed values; this way, a potential extreme effect of one high value was avoided (Byrnes and others 2014). For the threshold approach, the same maximum value was used and the threshold based indices were calculated as a sum of the number of functions surpassing a certain percentage of the maximum observed value. We used multiple single thresholds; 25, 50, and 75%.
The two approaches highlight slightly different perspectives of multifunctioning; the averaging approach helps to determine whether the average level of multiple functions changes with biodiversity. However, it does not distinguish whether functions are being contributed at similar levels, or if one function is very high and another function is low (that is, functions can compensate for each other: Byrnes and others 2014). Functional compensation was unlikely in this study because there were no strong negative correlations between the functions (Appendix 2). The threshold approach, on the other hand, shows the exact number of functions that exceed a certain threshold in relation to levels of biodiversity; therefore, the included functions cannot compensate for each other. We used these two multifunctionality approaches to explore relationships between ecosystem functionality and biodiversity (including species richness, abundance, and functional richness) along stressor gradients.

Statistical Analysis
We expected that there would be break points in the relationships between the ecosystem functions and the light availability, that is, there would be nonlinear responses to changes in light especially for the sediment primary production. We used regression trees (RPART-package in R; Therneau and others 2014) to objectively determine the potential break points for the four included functions (SOC, GPP, Chl a/phaeo-ratio, and CD) associated with both stressors (mPAR and porewater ammonium concentration). All analyses were conducted on the entire dataset that included all replicates from the nutrient treatments. Regression trees revealed that light (lmol photons m -2 s -1 of PAR) was forming the first splits for all four functions, and no first splits based on porewater nutrient concentrations were observed (all except one of the further splits in the regression trees were also formed due to light). The break points based on light were at 150 lmol photons m -2 s -1 of PAR for SOC, 344 for CD, 350 for Chl a/ phaeo-ratio, and 424 for GPP. To ensure that the responses were truly nonlinear we also performed linear regressions, but no significant relationships were identified. Hence, based on the regressiontree breakpoints, the data were split into three light groups: clear (> 420 lmol photons m -2 s -1 of mPAR; n = 63 samples), turbid (< 340 lmol pho-tons m -2 s -1 of mPAR; n = 72 samples), and moderately turbid (340 lmol photons m -2 s -1 £ mPAR £ 420 340 lmol photons m -2 s -1 ; n = 70 samples).
In order to elucidate the effects of the two stressors on the macrofauna communities and the relationship between the macrofauna and multifunctioning, we performed a suite of univariate analyses. First, we ran ANOVAs on the univariate macrofauna diversity measures S, N, and F Ric , with light group (3 factor levels: clear, moderate, turbid), nutrient addition treatment (3 factor levels: control, medium, high), and their interaction as independent variables. These tests were performed in order to determine how the two stressors interactively affected macrofauna communities. Second, we ran multiple linear regressions with the multifunction index (MF index ) as the response variable. These regressions included light group (factor with 3 levels), nutrient addition treatment (factor with 3 levels), univariate measures of diversity (continuous variable S, N, and F Ric ; one at a time) and the variables' interactions as predictors. If any interaction was significant, the relationships were further analysed within the groups, and an ANCOVA-type of analysis was run to reveal significant trend differences (that is, differences in slopes and intercepts) between groups.
For the analyses with the multifunctionality thresholds as the response variables, the same approach was used but generalised linear models (glm) with Poisson error distribution and a log-link function were run. Poisson error distribution was used because the response variables could only be integers, that is, 1-4 functions. Type-III sum of squares were consistently used; hence, if the higher order interaction was not significant (p > 0.15), it was removed and the analysis run again. All statistical analyses were performed in R (RStudio Team 2019).

Site Environmental Characteristics
The sites encompassed a gradient of water column turbidity, with light availability at the seafloor ranging from 134 to 712 lmol m -2 s -1 of submerged daily average PAR (mPAR; Table 1). The sites also encompassed a range of sediment characteristics: the ranges of the sediment median grain size and mud content (% < 63 lm) were 84-349 lm and 0-25%, respectively, with the sediment organic matter content varying between 0.7 and 5%. On average OM and mud content were higher within the moderately turbid and the turbid light groups compared to the clear light group (Table 1). The sediment nutrient enrichment successfully increased the porewater ammonium concentrations. The increase relative to the control porewater concentrations at the sediment surface (0-2 cm) was on average 40 times higher within the medium treatments and 250 times higher within the high treatments, whereas, deeper in the sediment (5-7 cm), the increase in medium treatment was on average 50 times higher and in high treatment 170 times higher compared to the control concentrations (Table 1; Appendix 3). On average, the natural porewater concentrations of ammonium (control PW, Table 1) were also higher within the moderately turbid and turbid sites compared to the clear sites.

Fauna Community Response to the Stressors
Overall, 122 species were observed across all sites, with an average of 22 species occurring at each site within the control plots (see Appendix 4 for a list of fauna) and a site-average abundance of 93 (range: 11-319) individuals per 0.01 m 2 in the control plots. The most numerous species were the bivalve Austrovenus stutchburyi, the spionid polychaetes Aonides trifida and Prionospio aucklandica, while the species occurring in the most samples were the bivalves A. stutchburyi (occurrence in 100% of the samples) and Macomona liliana (96% of the samples), and the capitellid polychaete Heteromastus filiformis (84% of the samples). Species richness was significantly affected by both the nutrient treatment and light group, but there was no significant interaction between these two stressors (p > 0.15, ANOVA; Appendix 5). The average number of species was 18 in clear sites compared to the moderately turbid and turbid sites with 21 and 22 species on average, respectively ( Figure 2; Table 1). All light groups had however similar ranges for species richness (Table 1). The average number of species in the high nutrient treatment was 18, compared to the control and medium treatments with 22 species on average in each (Figure 2). The same effects of treatment and light group were found for abundance and functional richness (Appendix 5).

MF index
Multiple linear regression showed no significant effect of the nutrient addition on MF index (Table 2).

Stressors Increase the Impacts of Coastal Macrofauna…
Declining species richness resulted in a decrease in MF index across the sites, but importantly this relationship was dependent on light group (significant light group 9 species richness interaction in linear mixed models; Table 2). The MF index decreased with decreasing number of species in turbid and moderately turbid sites, but there was no significant relationship within the clear light group (Figure 3A; Table 2). The biodiversity-ecosystem functioning relationships were thus nonlinear in response to turbidity. The slopes and intercepts of the regression lines within the turbid and moderately turbid groups did not significantly differ from each other. Similar results were found when using abundance or functional richness instead of the number of species as the predictor representing the benthic macrofaunal community (see Table 3, Appendix 6 and 7).

Multifunctionality Threshold Index
GLMs revealed that as species richness decreased, the number of functions exceeding the 25 and 50% thresholds also decreased; however, again this was dependent on the light group (significant light group 9 species richness interaction; Table 2). Declining species richness resulted in decrease in the number of functions exceeding the 25% threshold only in the turbid sites ( Figure 3B). In contrast, the declining species richness caused a decrease in the number of functions exceeding the 50% threshold in both the moderately turbid and turbid sites ( Figure 3C). In the clear sites, species richness was not significantly related to changes in the number of functions above the thresholds. A declining species richness resulted in a decrease in the number of functions exceeding the 75% threshold regardless of the light group (Table 2; Figure 3D). Similar trends were found when relationships between the thresholds and abundance and functional richness were examined (see Table 3, Appendix 6 and 7).

DISCUSSION
The need to examine the impacts of multiple stressors within natural variability was already highlighted two decades ago by Breitburg and  . This study incorporated a multistressor experiment into substantial real-world heterogeneity (wide range of latitudes and environmental contexts) to investigate the combined effects of two common estuarine stressors; turbidity and nutrient enrichment. The overall trends of the biodiversity-ecosystem multifunctioning relationships were nonlinear in response to turbidity, but interestingly no relationships were observed with the nutrient additions. Only a few significant relationships between multifunctionality and biodiversity were found within the clear sites, whereas in the moderately turbid and turbid sites the multifunctionality increased with increasing biodiversity, regardless of the index used (Table 3). This would suggest that a decrease in number of species, abundance or functional richness may result in a decline in the ecosystem functioning in moderately turbid and turbid systems. Whereas, in the clear systems the multifunctioning remained more constant regardless of decreasing biodiversity, which suggests that these clear systems are more resilient towards a loss of biodiversity compared to turbid systems. The overall ecosystem functioning in more oligotrophic conditions not stressed by turbidity is potentially underpinned by more pathways between species, environmental properties and functions compared to the more turbid and nutrient rich conditions. Figure 2. Species richness as a function of A the light groups and B the nutrient addition treatment group. Lower case letters above the boxes indicate significant differences among groups (p < 0.05). Outliers are defined as > 1.5 times the interquartile range beyond either end of the box. Thrush and others (2021) showed through ecosystem interaction network analysis that within clear systems, the connectivity was high between the different ecosystem components, abiotic and biotic units with oxygen and nitrogen fluxes, and denitrification, compared to simpler networks in turbid systems.
The way an ecosystem loses species due to stress is not random, and community responses to stress may vary depending on the starting point and which community is remaining (Zavaleta and Hulvey 2004;Srivastava and Vellend 2005;Bracken and others 2008;Siwicka and others 2020). A more diverse community is thought to be more resilient against disturbance, due to the more possible responses the collection of species can have (for example, Mori and others 2013). Marine benthic communities typically contain only a few very numerous species, while most species are only represented by a few individuals (Ellingsen and others 2007), and these rare species are likely more vulnerable to stress (Walker 1992;Gray 1997;Jackson 2001). In this study, an examination of the ratio between common and rare species among sites showed that on average, the alpha diversity was higher in the turbid sites, while the clear sites had a higher beta diversity. The most species rich samples (> 90th percentile) within the clear sites contained 28% rare species, compared to the species rich samples in turbid containing only 12% rare species (rarity defined based on species occurrence < 25th percentile and abundance < 25th percentile within the control samples of each light group). Within the clear sites, we also lost a larger proportion of rare species (43%) when going from high-diversity samples to lowdiversity samples compared to the turbid sites (22%). In clear light conditions as species richness was reduced, a lower number of common species were lost and since the common species likely contribute most to functioning, this might partially explain why multifunctionality remained constant with the biodiversity loss (Table 2). No interaction between the three predictors was significant (p > 0.15) therefore not reported. *pseudo-R 2 is from Nagelkerke (Cragg and Uhler); values are between 0-1, that is, similar interpretation to the R 2 from linear models.
Common and numerous species are typically used in biodiversity-ecosystem functioning studies, and they are often the most important for functioning (Ellingsen and others 2007). Common species have even been reported to compensate for losses in rare species, at least on short time scales and with specific observed functions (Smith and Knapp 2003). The results of this study would suggest that the more common species contributed to the measured functions, hence the decline in functioning at the turbid sites where more common species were lost. One should however remember that rare species, due to their unique set of attributes (Ellingsen and others 2007;Mouillot and others 2013), might contribute to ecosystem resilience in other ways that are not measured here, for example turnover over longer time scales and when other functions are considered (Loreau and others 2001;Hector and Bagchi 2007;Gamfeldt and others 2008;Siwicka and others 2021). The maintenance of function with species loss in the clear sites could also signal higher functional redundancy and therefore an ability to lose a level of biodiversity without impacts on function.
To understand how multiple stressors affect biodiversity and multiple ecosystem functions, different methods, such as monitoring, modelling, and experimentation in a wide range of environments are going to be needed. Here, we incorporated a large-scale field experiment with natural heterogeneity to examine biodiversity-ecosystem functioning responses to two common stressors, and we observed a potential indirect interaction between the stressors that can lead to a decline of ecosystem function. Our study found that nutrient addition negatively impacted species richness, and whilst this did not directly interact with turbidity, we found that turbid environments were more vulnerable to functional losses with declines in species richness. This would imply that estuaries already stressed from high water column turbidity are likely more vulnerable to increased levels of nutrients due to the possible negative effects on benthic macrofauna communities, which in turn have a negative effect on functioning. These results highlight an overarching snowballing effect that cumulative effects can induce and this kind of mechanistic understanding is important for informing predictions of abrupt ecological shifts. Figure 3. Examination of the relationships between A the average multifunction index (MF index ) and number of species within the light groups, B the number of functions exceeding the 25% threshold, C the 50% threshold within the light groups, and D the number of functions exceeding the 75% threshold and the number of species. Only significant relationships illustrated. There was no significant relationship within the clear light group for either approach, but similar trends of increasing multifunctioning with increasing species richness were found in the moderately turbid and turbid light groups (Table 2). For the 75% threshold, only number of species, independently of light group, was significant.

AC KNOWLEDGEMENT S
We thank the large consortium of colleagues around New Zealand for their help in collecting all this extensive field data set. The project was funded through the New Zealand National Science Challenge Sustainable Seas, Dynamic Seas, Tipping Points project (CO1X1515).

FUNDING
Open Access funding enabled and organized by CAUL and its Member Institutions.

De clarations
Conflict of interest The authors declare that they have no conflict of interest.

OP E N A C C ES S
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/. Table 3. Summary of the Significant Effects (p < 0.05) of Species Richness (S), Abundance (N), and Functional Richness (F Ric ) on the MF index and the Multifunctionality Thresholds and Their Directions Indicated (see Table 2, Appendix 6 and 7 for details)