Cumulative stressors reduce the self-regulating capacity of coastal ecosystems.

Marine ecosystems are prone to tipping points, particularly in coastal zones where dramatic changes are associated with interactions between cumulative stressors (e.g., shellfish harvesting, eutrophication and sediment inputs) and ecosystem functions. A common feature of many degraded estuaries is elevated turbidity that reduces incident light to the seafloor, resulting from multiple factors including changes in sediment loading, sea-level rise and increased water column algal biomass. To determine whether cumulative effects of elevated turbidity may result in marked changes in the interactions between ecosystem components driving nutrient processing, we conducted a large-scale experiment manipulating sediment nitrogen concentrations in 15 estuaries across a national-scale gradient in incident light at the seafloor. We identified a threshold in incident light that was related to distinct changes in the ecosystem interaction networks (EIN) that drive nutrient processing. Above this threshold, network connectivity was high with clear mechanistic links to denitrification and the role of large shellfish in nitrogen processing. The EIN analyses revealed interacting stressors resulting in a decoupling of ecosystem processes in turbid estuaries with a lower capacity to denitrify and process nitrogen. This suggests that, as turbidity increases with sediment load, coastal areas can be more vulnerable to eutrophication. The identified interactions between light, nutrient processing and the abundance of large shellfish emphasizes the importance of actions that seek to manage multiple stressors and conserve or enhance shellfish abundance, rather than actions focusing on limiting a single stressor.

Introduction

Humans have settled along coastal margins and estuaries throughout history because they provide vital food resources and transport routes. However, major past and ongoing changes to coastal ecosystems have resulted from the development of ports, cities, agriculture, forestry, fishing, and industrialization (Nichols et al. 1986, Cooper and Brush 1993). Despite the small proportion of the marine environment classified as coastal or estuarine (~4% land area and 11% ocean [Barbier 2017]), these ecosystems are crucially important for global carbon and nitrogen cycling and other ecosystem services yet are among the most impacted. With the increase in global human populations around coastal margins, biodiversity loss and climate change impacts have escalated, predisposing these systems to risk of abrupt nonlinear shifts in ecosystem functioning (Lotze et al. 2006, Barbier et al. 2008, Cloern et al. 2016). Tipping points move a system between functionally different states, with reversal difficult due to feedbacks that lock a system in a degraded state (Nyström et al. 2012). Ecological tipping points can occur through large changes in external factors (climate change, earthquakes, or resource extraction) or, more challengingly, through subtle and gradual changes that influence interactions between ecosystem components (Carpenter et al. 2001, Scheffer et al. 2001, Hastings et al. 2018).

Predicting abrupt change in real‐world ecosystems remains a major challenge for ecology and, while improving the quality of long‐term data is essential, it alone is insufficient to allow us to develop the understanding that will allow management actions that prevent rapid transitions (Ratajczak et al. 2018, Hewitt and Thrush 2019). Identifying the existence, location, and cause of threshold changes is most clearly achieved by experimentation. Conducting experimental manipulations on chronic and accumulating stressors is difficult, requiring targeted studies merging ecological theory with correlative and manipulative studies. Such studies may prove highly effective in driving insight into appropriate coastal management strategies. These strategies are required as, despite the need to address cumulative effects in management being well‐recognized and a number of frameworks existing (Crain et al. 2008, Ban et al. 2010, HELCOM 2017, Stelzenmüller et al. 2018), the use of additive effects models to define environmental risks is widespread (Furlan et al. 2019). Interactions between stressors that change across thresholds is a particularly important topic for coastal ecosystems containing high proportions of soft sediment habitats, within which strong interactions generally exist between the resident species and the physical and chemical environment.

In estuarine and coastal environments, shellfish beds influence many important functions that maintain ecosystem health. Bivalves that protrude from the sediment affect boundary flows and provide refugia and settlement surfaces for other species, thus enhancing biodiversity (Hewitt et al. 2005). Suspension feeding bivalves play crucial roles in actively filtering (and clearing) the water column and depositing organic rich biodeposits on the sediment surface (Norkko et al. 2006). Bivalves that live in the sediment move and mix sediment particles, with species at the sediment–water interface (e.g., Austrovenus stutchburyi) indirectly changing sediment topography via bulldozing and biodeposit production (Woodin et al. 2016). Subsurface species (such as Macomona liliana) influence nutrient cycling by generating pore water pressure waves that produce rapid redox oscillations deep in the sediment (Volkenborn et al. 2012). Thus, at sufficient abundance and biomass, shellfish can profoundly alter how coastal ecosystems function through changes to water quality, sediment biogeochemistry and, in particular, the remineralization of organic matter, and the flux of nutrients and oxygen across the sediment water interface. As a major fraction of the increasing load of nitrogen entering the coastal zone is denitrified in coastal sediments (Middelburg et al. 1996, Galloway et al. 2008), it is likely that shellfish play a critical role in regulating ecosystem responses to eutrophication.

However, shellfish beds are critically degraded in many coastal habitats (Airoldi and Beck 2007, Beck et al. 2011), not only impacted by common stressors (eutrophication, sedimentation, and habitat disturbance) but also by human harvesting. The snowballing effects of the functional extinction of shellfish beds coupled with elevated nutrient and sediment delivery to coasts may lead to rapid changes disproportional to nutrient loading, predisposing the ecosystem to a functional tipping point. As a result, management actions need to be based on an understanding of the complexities, including nonlinear interactions, of multiple stressors on overall ecosystem functioning.

There are a wide variety of model‐based approaches that have been developed to identify tipping points and thresholds and inform management and policy (Samhouri et al. 2010, Kelly et al. 2015), (Samhouri et al. 2017). These are typically applied to detect thresholds after the event (Carpenter et al. 2001, Clements and Ozgul 2018). One strategy for unravelling complexities in real‐world ecosystems is to develop ecosystem interaction networks (EIN) (Thrush et al. 2014). EINs are models based on connections between ecosystem components, although generally only a subset of components that represent a full ecosystem are included. EINs allow exploration of feedbacks, direct and indirect effects and modifiers, all of which are important factors for understanding the potential for cascading effects and self‐regulation. Thus, EINs allow us to shift our focus from recognizing simple cause–effect relationships to understanding interactions across critical ecosystem components and can be used to identify key connections upon which a network is dependent. These key connections may be state dependent, i.e., may change as a threshold is crossed, thus, combined with approaches that specifically identify threshold responses between variables, EINs can be a powerful tool to demonstrate environmental risks.

In this study, we defined a hypothesized EIN for the functioning of shallow coastal areas involving shellfish, denitrification, ammonium and oxygen fluxes between the sediment and the water column, benthic primary production, sediment grain size, nitrogen, light, and temperature (Fig. 1 and Appendix S1). The connections between, and the relative importance of, components of the EIN were predicted to change with chronically elevated turbidity, particularly components related to the ability to process nitrogen with increasing nitrogen. This hypothesis was tested with data from a long‐term large‐scale field experiment. This experiment allowed us to assess how the effect of nitrogen additions (experimental manipulation), site‐dependent variations in turbidity and the cumulative effects of the two stressors influenced the architecture of the EIN. Pore‐water nutrient concentrations were manipulated at sites arrayed along a turbidity gradient (as a proxy of suspended sediment concentrations driven by terrestrial soil erosion). The experiment spanned 10° of latitude at 24 sites in 15 estuaries, allowing us to generalize cumulative effects and determine real‐world thresholds at a national scale. Specifically, we predicted that there exists a threshold in the amount of light available for photosynthesis at the seafloor (associated with suspended sediment) below which (1) the ecosystem interaction network would be relatively simple, (2) the importance of shellfish in the processing of nutrients would be reduced, and (3) the ecosystem would be less able to cope with increased nutrients. That is, rather than nitrogen and turbidity being additive stressors, there would be an interaction between them that included a threshold.

Fig. 1

Hypothesized interaction network linking shellfish (Austrovenus stutchburyi and Macomona liliana) to sediment properties and ecosystem functions (see Appendix S1: Table S1 for details and references of the connections). To differentiate the types of variable, teal background indicates external drivers; red indicates shellfish; purple indicates sediment (mud, particles < 63μm in sediment; gravel, particles >2 mm in sediment [generally shell hash]; OM, organic matter); blue indicates benthic oxygen consumption (BOC), green indicates plants (benthic chl a, standing stock of microphytobenthos); yellow indicates nitrogen (den, denitrification (N2); NH4 flux, flux across the sediment–water interface.

Materials and Methods

Study design

Using local knowledge, we identified 24 sites in 15 estuaries across New Zealand (Fig. 2) that encompassed a wide variation in water turbidity. Utilizing a latitudinal gradient of 10 degrees we designed the experiment to obtain a turbidity gradient that was not spatially confounded. Data collected in the experiment (Data collection) confirmed that latitude was not a strong driver of turbidity or ambient nutrient conditions (control pore water ammonium concentrations) at our sites, with all Pearson’s (and Spearman's) correlations coefficients between these three variables being <0.25).

Fig. 2

New Zealand with the locations of the 15 study estuaries marked (some estuaries had multiple sites with a total of 24 experimental sites). WRG, Whangarei Harbour; WGT, Whangateau Harbour; MAH, Mahurangi Harbour; WHI, Whitianga Harbour; RAG, Raglan Harbour; TAU, Tauranga Harbour; DEL, Delaware Inlet; WMA, Waimea Inlet; AVO, Avon‐Heathcote Estuary; AKA, Akaroa Harbour; BLU, Blueskin Bay; WKW, Waikawa Estuary; NEW, New River Estuary; JAC, Jacobs Creek Estuary.

Sites were situated at about mid‐tide level, on sandy permeable sediments and in areas occupied by the functionally important bivalves Macomona liliana and Austrovenus stutchburyi at their natural densities (see Results). At each site we established three 9‐m2 plots for each of three treatments, high fertilizer 600 g N/m2, medium fertilizer 150 g N/m2, and disturbance controls. We used Nutricote (Chisso‐Asahi Fertilizer Co., Ltd, Tokyo, Japan) slow release fertilizer (40‐0‐0 N:P:K) injected uniformly into the sediment at a depth of 15 cm to elevate porewater nitrogen concentrations for extended periods (Douglas et al. 2016, Douglas et al. 2018). In the context of our experiment, we are not using the nutrient treatments to test for treatment differences as we would in a classical hypothesis testing experiment, rather we are using the treatments to investigate the ecosystem’s nutrient response in the context of the EIN typology. The experimental plots throughout the country were set up March–April 2017.

Data collection

At each site for the duration of the experiment, we deployed water column photosynthetically active radiation (PAR) sensors recording every 10 minutes (fixed vertically 10 cm above sediment surface; Odyssey, Dataflow Systems, Christchurch, New Zealand) and sediment temperature sensors logging temperature every 30 minutes (buried at depths of 3 and 7 cm; i‐Buttons, Maxim Integrated Products, San Jose, California, USA) to establish records of environmental change over the course of the experiment. The plots were occasionally checked to service loggers but left undisturbed for the next 7 months. After downloading, the PAR readings during submersion were summarized into daily averages, maximums and minimums. Averages of these were created for each site (ADAL, average daily average; ADML, average daily maximum; and ADLL, average daily minimum) for use within the statistical analyses.

The experimental plots were sampled for macrofauna, sediment characteristics, and biogeochemical fluxes in October–November 2017 (springtime in New Zealand). We used multiple teams around the country to optimize sampling relative to the timing of tides. Five sediment cores (2.6 cm diameter, 2 cm depth) were collected and pooled from each plot for analysis of sediment chlorophyll a, grain size, and organic content, and stored frozen (−20°C). Sediment grain size was measured with a Malvern Mastersizer‐3000 (Malvern, UK) (particle size range 0.01 to 3,500 µm) after removal of organic matter with 10% hydrogen peroxide (Singer et al. 1988). Chlorophyll a and phaeopigments were extracted from sediment in 90% buffered acetone and measured before and after acidification using a Turner Designs (San Jose, California, USA) 10‐AU fluorimeter (Arar and Collins 1997). Sediment organic matter content was determined by weight loss on ignition (60°C for 24 h or until stable mass + combustion for 4 h at 550°C; Dean 1974).

Two macrofauna cores (13 cm diameter, 15 cm deep) were collected from each plot, sieved (500‐µm mesh), and preserved in 70% isopropyl alcohol. Macrofauna were stained with Rose Bengal, sorted, and identified. Macomona and Austrovenus were measured (longest shell dimension) and the numbers of individuals sized >20 mm were totaled, for each species, across the two replicate cores for each plot. This size of bivalves was selected to represent the large adults that previous studies have shown to play important roles in sediment functioning (Thrush et al. 2006, Jones et al. 2011, Thrush et al. 2014, Woodin et al. 2016).

Benthic flux chambers were used to determine biogeochemical fluxes. We had multiple teams operating around the country to minimize variation in sampling dates between locations (26 October–27 November 2017) while ensuring that chambers were deployed on sunny days with mid‐day high tides. Hobo temperature loggers (Onset HOBO Pendant Temperature/Light Data Logger, Bourne, Massachusetts, USA) were deployed in chambers in case temperature was required for the EINs. Biochemical measures made were net sediment–water O2, N2, and NH4 fluxes. NOx and P fluxes were also measured but as a large proportion of the data was close to detection limits, these variables were not included in our analyses. In each plot, we deployed two chamber bases (50 × 50 × 15 cm height) pressed 5 cm into the sediment during low tide. On the incoming tide (water depth ~50 cm), acrylic dooms sealed ~40 L of ambient seawater over the sediments. Opaque shade cloth prevented light entering one chamber from each plot. Chambers were incubated for ~4 h over a midday high tide. Water samples (1 × 60 mL syringe for solute, and 2 × 60 mL airtight syringes for gas concentrations) were withdrawn from the chambers through sampling ports at the beginning and end of the incubation period. Dissolved O2 concentrations in the water samples were measured using an optical probe (ProODO YSI, Yellow Springs, Ohio, USA). Samples were then filtered through a 0.45‐µm Whatman GF/C glass fiber filter and stored frozen (−20°C) prior to analysis of ‐N Water from airtight syringes were preserved in zinc chloride and stored in gas tight exetainers (Labco, Ceredigion, UK) at 4°C until analysis of N2 concentration. Solute and gas fluxes were calculated as (C end− C initial × V)/A × T, where C is nutrient or oxygen concentration (µmol·L−1·L−1), V is the volume of seawater inside the chamber (L), A is the area of sediment enclosed by the chamber (m2), and T is the elapsed time between initial and final samplings (h).

Water samples collected from the chambers were analyzed for with a Lachat QuickChem 8000 Series FIA+ (Zellweger Analytics, Milwaukee, Wisconsin, USA) using standard operating procedures for flow injection analysis (NH4 +‐N detection limit of 0.07 μmol/L). N2 gas fluxes were analyzed on a quadrupole membrane inlet mass spectrometer (MIMS; Bay Instruments, Cambridge, Maryland, USA) using the N2/Ar method (analytical precision <0.03%; [Kana et al. 1994]). Loss of samples in transit meant that data for MIMS analysis were lost from the JAC, NEW, WKW, and BLU sites.

EIN analyses

While we had a priori predicted that there would be a primary threshold associated with light, we tested for nonlinear relationships between our ecosystem components and both light and nutrients. Regression trees were then used to objectively determine the position of any break points associated with our predictor variables (submerged daily average PAR reading [averaged across the experimental duration; hereafter called ADAL] and nutrient treatment) for three of the ecosystem components thought to be most likely to be directly affected (denitrification, NH4 efflux, and benthic oxygen consumption). Regression trees explain variation in a single response variable by repeatedly splitting the data into two more homogeneous groups, using the best explanatory variable in each case, and thus the break point delineates a (significant) change in the relationship between the stressor and the ecosystem function. We also tested to ensure that the regression tree gave better results than a linear regression (i.e., we truly did have a nonlinear response). Linear regressions detected no significant relationships and explained <10% of the variance. However, the regression tree results were all significant and varied between explaining 20% to 33% of the variance, suggesting that relationships were nonlinear and the breakpoints detected were useful.

Regression tree analyses were conducted using the RPART package (Therneau et al. 2014). Tree growth was constrained to have a minimum of 20 observations in a node (group) before attempting a split; the split had to increase the fit (represented by the R 2) by ≥0.03 and each terminal node (final most homogeneous group) had to contain at least 10 data points (i.e., two sites). Cross validation and tree pruning were not used.

ADAL was the most important driver of variability for all three of the variables, forming the first tree split (Table 1). Regression tree analyses were then conducted with only ADAL as a predictor variable. These analyses demonstrated the majority of break points occurred at or between 350 and 420 μmol·L−1·m−2·s−1 PAR (Fig. 3). Once we knew that the regression tree analysis not only suggested a nonlinear response but that the suggested break points for PAR were similar across the selected ecosystem functions, we could then go on to analyze whether there was a difference between the networks associated with the suggested break points. Therefore, our data were broken into three groups: Clear sites with >420 μmol·L−1·m·2·S−1 PAR, turbid sites with <350 μmol·L−1·m−2·S−1 PAR; and an intermediate set with too few data points to warrant further analysis.

Table 1

Results of regression trees used to objectively determine the position of break points for three variables (denitrification, ammonium efflux, and benthic oxygen consumption) with two explanatory variables, average daily light when the sites were inundated (ADAL) and nutrient treatment.

Variable	ADAL (μmol·L−1·m−2·S−1 PAR)	Nutrient addition
Denitrification	350 (1)	>control (4)
NH4 flux	417 (1)	<high (2)
Benthic oxygen consumption	340 (1)

The values or treatments at which a split first occurred are indicated along with the split level (in parentheses).

Fig. 3

Regression trees used to objectively determine the position of any break points for three variables (denitrification, ammonium efflux, and benthic oxygen consumption) related to the average daily average photosynthetically active radiation (PAR) when the sites were inundated (ADAL).

EIN derivation using Structural Equation Models (SEM)

SEM (Grace et al. 2010, Kline 2011) were used for our EIN as this approach starts with the definition of the potential interaction network, built using expert opinion and literature (Fig. 1 and Appendix S1: Table S1). Our network encompassed the main linkages between sediments, nutrients, bivalves, and light. SEM analysis was conducted on the clear and turbid groups separately; the intermediate group was not analyzed as it contained few data points. Information from the three replicate plots of each treatment at each site was included as plot‐specific data on all variables (with the exception of ADAL) was available. SEMs were developed on standardized data using M‐Plus software (Muthén and Muthén 2007). Data transformations (log) to improve the linearity of responses, or normality of errors, were needed for benthic oxygen consumption and average daily maximum PAR. The hypothesized EIN was pruned for the clear and turbid data separately by removal of direct pathways based on parameter tests of significance (attempting to remove pathways with P > 0.10), while maintaining non‐significance tests for overall SEM goodness of fit (χ2) and the root mean square error of approximation, increasing values of the comparative fit index (above 0.95; Vile et al. 2006) and decreasing values of Akaike’s Information Criterion (Table 2). We compared the two pruned EINs based on differences in ecosystem components retained, network typology, and the directionality and strength of relationships between ecosystem components.

Table 2

SEM results for the best model obtained for the clear and turbid data sets.

Statistic	Clear	Turbid
Akaike information criterion (AIC)	1,162.1	883.5
Bayesian information criterion (BIC)	1,251.4	933.2
Sample‐size adjusted BIC (n* = (n + 2)/24)	1,119.2	863.9
χ2	31.9	13.72
P	0.0598	0.6867
RMSEA	0.09	<0.01
P	0.149	0.824
CFI	0.951	0.999

The χ2 is for goodness of fit (P < 0.05 indicates model significantly different from the empirical data); RMSEA, root mean square error of approximation (close to zero is good fit and P tests for difference from zero); CFI, comparative fit index (ranges 0–1, with 1 being a perfect fit).

We also determined whether there was real separation between the two pruned EINs by testing the significance and goodness of fit of the turbid EIN structure on the clear data and vice versa. Imposing the turbid EIN structure on the clear data demonstrated that the empirical data could not support the model. Imposing the clear EIN structure on the turbid data similarly demonstrated that a lack of fit between the empirical data and the model.

Results

The experiment encompassed a wide range in each of the parameters used to develop the SEMs, although the interquartile ranges were more constrained (Fig. 4). Model pruning resulted in an SEM for the turbid data with a very good comparative fit index and no pathways having P values > 0.10 (Tables 2, 3). However, for the clear data, achieving a reasonable fit required retaining a number of pathways of low significance (Table 4). Removal of these pathways, including the variables predicted to affect NH4 flux, resulted in χ2 values and RMSEA with P values ≪0.05 and CFIs below 0.9. Thus, these pathways, while weak, are important.

Fig. 4

Box plots showing the median (mid line), interquartile range (box edges), and range (whiskers) for sediment, shellfish, and ecosystem function variables used in the structural equation model (SEM).

Table 3

Turbid SEM parameter estimates and two‐tailed tests for whether the estimate differs from zero.

Predicted and explanatory effect	Estimate	P
Sediment–water flux NH4
Benthic oxygen consumption	+0.03	0.0460
Austrovenus	−0.03	0.0520
Benthic oxygen consumption
Chlorophyll a	+0.40	0.0010
Chlorophyll a
Austrovenus	+0.62	<0.0001
Denitrification	+0.19	0.1000
Austrovenus
Mud	+0.42	0.0150
Macomona
Mud	−0.27	<0.0001
Mud
Organic matter	+0.71	<0.0001
Chlorophyll a	−0.50	0.0010

Table 4

Clear SEM parameter estimates and two‐tailed tests for whether the estimate differs from zero.

Predicted and explanatory effect	Estimate	P
Denitrification
Benthic oxygen consumption	+0.46	0.0300
Average daily maximum PAR	−0.42	0.2400
Organic matter	+0.67	0.2340
Sediment–water flux NH4
Austrovenus	+7.86	0.7860
Benthic oxygen consumption	−13.49	0.7800
Macomona	−2.13	0.7640
Mud	+7.23	0.7760
Benthic oxygen consumption
Austrovenus	+1.07	0.0210
Sediment–water flux NH4	+2.88	0.0380
Chlorophyll a	−1.26	0.0860
Average daily maximum PAR	+1.64	0.0730
Macomona	+0.55	0.1760
Experimental N addition	−0.55	0.1490
Average daily maximum PAR
Chlorophyll a	−0.07	0.5830
Chlorophyll a
Austrovenus	−1.21	0.3160
Denitrification	−0.75	0.3720
Average daily maximum PAR	+0.78	0.5100
Mud	−1.75	0.4670
Austrovenus
Chlorophyll a	+0.38	0.1960
Macomona	−0.13	0.1980
Macomona
Mud	+1.59	<0.0001
Organic matter	−0.93	0.0130
Mud
Chlorophyll a	+0.65	<0.0001
Organic matter
Mud	+0.74	<0.0001

PAR, photosynthetically active radiation.

The EINs derived from the experimental data demonstrate distinct differences in ecosystem structure (the components remaining as drivers) and connectivity between clear (>420 μmol·L−1·m−2·s−1 PAR) and turbid (<350 μmol·L−1·m−2·s−1 PAR) water sites (Fig. 5). Most importantly, the EIN derived from the clear sites is more complex, containing 22 connections (including two‐way links), whereas the turbid SEM contained only 9. In both cases, there are a similar ratio of positive to negative relationships (clear: 11:11, turbid 4:5, +:−).

Fig. 5

Ecosystem interaction network from (A) clear: >420 μmol·L−1·m−2·s−1 PAR and (B) turbid sites: <350 μmol·L−1·m−2·s−1 PAR. Teal boxes relate to external factors (nutrient addition is experimental manipulation; Ave. daily max PAR is the effect of turbidity on available photosynthetically active radiation). Purple boxes denote sediment characteristics (Chl a represents the standing stock of microphytobenthos; mud, mud in sediment; OM, organic matter in sediment). The blue box relates to oxygen (BOC, benthic oxygen consumption). Yellow boxes relate to nitrogen processing (Denitrification (N2); NH4 flux, flux across the sediment–water interface. Pink boxes relate to the abundance of large shellfish (individuals >20 mm shell length). Red arrows indicate positive relationships and black negative. Arrow thickness indicates relationship strength.

The experimental nutrient addition and the role of light were decoupled from any components in the EIN across the turbid sites, which together with the lack of any predictive connections in the model to denitrification, highlights reduced capacity of the ecosystem to process nutrient additions. For the clear data, however, denitrification was driven directly by benthic oxygen consumption and organic matter (positive) and average daily maximum PAR (negative; Table 3). Denitrification was reduced as nitrogen increased through the indirect effect on reduced benthic oxygen consumption. The strongest effect on benthic oxygen consumption was by sediment‐water flux NH4 (highest standardized parameter estimate, Table 4), followed by average daily maximum PAR, chlorophyll, and Austrovenus.

The two shellfish (Macomona and Austrovenus) were more involved in the EIN at the clear sites and directly affected benthic oxygen consumption. Conversely, in turbid sites, only chlorophyll a, reflecting microphytobenthos standing stock, affected oxygen consumption (probably due to photosynthesis), and no components affected denitrification. Interestingly, we did not find an effect of nutrient addition on chlorophyll a in either clear or turbid water.

Discussion

We were able to identify functionally important changes in coastal soft sediments due to a threshold in elevated turbidity interacting with shellfish densities. Contrasts between the architecture of the two EINs highlight a reduced capacity of the sediments to process increased nutrients with elevated turbidity, which implies a snowballing effect where the ecosystem effects of nitrogen loading can be amplified by the reduced capacity of the sediments to process nitrogen. This result emphasizes the importance of not only managing loadings of multiple contaminants but also the ecosystem’s responses as observed through the interactions of key ecosystem components. Nitrogen loading and sediment runoff are coastal and estuarine stressors of global significance (Thrush et al. 2004, Diaz and Rosenberg 2008). Our results highlight how these stressors can interact to generate cascading effects in the loss of ecosystem function.

Our finding that incident light to the seafloor is a critical threshold variable affecting system functionality, particularly N cycling, supports previous studies emphasizing the importance of both the activity and standing stock of microphytobenthos on intertidal flats (Pratt et al. 2014, O'Meara et al. 2017, Douglas et al. 2018). The corroborative evidence of independent studies supports both the identification of thresholds and the importance of managing systems to limit the risk of crossing these thresholds. Previous research in New Zealand estuaries revealed shifts in network architecture across a benthic chlorophyll a threshold of 11.6 µg/g (Thrush et al. 2012), with a similar threshold detected in a field experiment that shaded the sediment–water interface (Thrush et al. 2014). Limiting the light available to benthic microphytes when intertidal flats are inundated affects the rates of organic matter and nutrient processing, as well as the primary source of productivity (microphytobenthos) that underpins many coastal food webs and ecosystem services (Hope et al. 2019). There has been a long history of using thresholds to inform environmental management, for example using light thresholds to determine the depth distribution or impact of elevated turbidity on seagrass (Duarte 1991) or using animal oxygen thresholds to determine species sensitivity to hypoxia (Vaquer‐Sunyer and Duarte 2008, 2011). These species‐based physiological thresholds are important but they are not the same as defining a threshold in the context of an EIN. EIN‐derived thresholds are developed in the context of cumulative effects and mechanisms of interactions between different ecosystem components. This provides a new tool for managers grappling to marry cumulative effects and ecosystem dynamics with a mechanistic understanding of the interactions between organisms, their environment and key stressors. The mechanistic information serves two critical functions in a management context: ensuring appropriate information on processes and rates are included in models relevant to management; and providing ecosystem natural history information that can be communicated across society to enhance ocean literacy. For example, these features come together in the context of managers coping with multiple stressor effects in estuaries. The EINs imply that national or regional standards for nitrogen loading are a blunt management tool in the face of variations in turbidity, but equally the EINs draw attention to the dynamics of ecosystem’s responses and their capacity for denitrification and processing nitrogen.

The shellfish in the EIN represent the control parameter that is likely to change the most slowly given the relative metabolic rates of microbes and macrofauna (see Appendix S1: Table S1), highlighting a potential for their functional extinction to induce significant change in ecosystem function and result in a loss of adaptive capacity. The nature of change in the EIN with increasing turbidity (that is, simplification and a dependence on a few strong connections) emphasizes the importance of interactions between physical and biological ecosystem components in changing ecosystem function. The application of complex system theory to the resilience of ecological systems has highlighted the importance of breaks in feedback loops in generating tipping points (Scheffer et al. 2001, Carpenter 2003, Scheffer 2009, Biggs et al. 2012, Selkoe et al. 2015). The differences we observed in EIN architecture indicates that soft sediments in turbid and clear water conditions operate in alternative states, albeit with feedbacks present in both systems. While this has been postulated from post hoc analysis of systems that have passed a tipping point (Nyström et al. 2012), our study shows, empirically, that there are changes in the wiring of the EIN associated with small changes across a turbidity gradient (Fig. 5). Such a result offers insight into how we can assist managers in developing more adaptive and case specific policy and management thresholds.

Our EIN analysis demonstrates how elevated turbidity (in our case associated with sediment runoff) decouples sediment nitrogen processes. As elevated turbidity reduces the capacity to denitrify and thus remove nitrogen, it becomes more likely that there will be a tipping point associated with eutrophication as the nitrogen builds up more quickly in the system. While we cannot “prove” this with the clarity and elegance of a theoretical model, or the hindsight of an extensive time series, we do provide empirical evidence, at a national scale, of important functional shifts in these complex biogeochemical systems. SEMs provide a way of analyzing empirical data to identify the architecture of networks of interactions, but they do have drawbacks, particularly related to pinning down the form of specific relationships, providing predictions and the amount of data required for systems with multiple connected components. They are also not temporally dynamic. In our case, New Zealand has a very temperate climate and does not experience the strong seasonality usually associated with continental landmasses. Thus, seasonality is not likely to be strong factor influencing our results. However, to really demonstrate this we would have to conduct our experiment at a number of sites at different times of the year to include both large differences in temperature and other likely seasonal drivers. Sediment temperature was not included in our final models (Fig. 5), because we sampled the experiment at a time of the year when there were only small differences between sites.

SEMs are not generally predictive models. Our SEMs were explanatory in nature, developed using covariance matrices and the maximum likelihood estimator. Partial least squares path modelling (PLSPM) can be used to create predictive SEMS in which case the composite model is developed using multiple regression (Hair et al. 2017). However, some dispute over the utility of PLSPM for prediction has been raised (Ronkko et al. 2016). Regardless, SEM provide an important and much needed link between empirical and theoretical research (Hastings 2016), explicitly testing whether theories about the degrees of connectivity and the importance of feedback loops matches the empirical data. In our experiment we focused on ecosystem components that are related to nutrient processing in the seafloor but ecosystems are multifunctional and other EINs could be constructed to address other functions (e.g., carbon sequestration, productivity, habitat formation) and these may have different tipping points.

The ability to identify networks of drivers of ecosystem function informs both increasing understanding of the degradation of estuaries, as well the potential for restoration and hysteresis to influence management interventions in degraded estuaries. Essentially management interventions in estuaries and coastal ecosystems can involve two actions: first, reducing stress loading through limit controls, spatial or temporal management, or marine reserve creation; and second actively restoring important species or habitats. Specifically, our results highlight that management objectives for turbid systems should consider a restoration focus to reduce multiple stressor effects, while in clear water estuaries a resilience focus could be developed to maintain the adaptive capacity of the ecosystem to process nitrogen and minimize sediment loading. Faced with nonlinear changes in ecosystems, making a decision around these actions requires information on where their specific system sits relative to a tipping point. Developing and testing appropriate EINs and collecting ecological and environmental information relative to changes in EIN architecture (e.g., Fig. 5) can provide managers with clues about whether a system is close to or has crossed a tipping point in ecosystem function. While this can be achieved when extensive data are available (Hunsicker et al. 2016), in many ecosystems relevant information is not available and in these circumstances experimental approaches can be used to gather useful data to test for changes in EINs. Moreover, the definition of restoration goals is critical for active restoration and these goals are often driven by the desire to restore ecosystem services. The EIN described in our study is focused on the role of infaunal shellfish in affecting sediment biogeochemistry, thus linking to services associated with productivity, biodiversity, water clarity and mitigating eutrophication. Active restoration of soft‐sediment habitats is still very much an emerging management practice but for seagrass, saltmarsh, and shellfish positive interactions and feedback processes are important for success (Grabowski et al. 2005, Suykerbuyk et al. 2012, Silliman et al. 2015), suggesting a role for SEMs in predicting restoration success.

Providing insight into the potential for non‐additive cumulative effects based on observed mechanisms of interactions and changes in feedback processes has important implication for ecosystem‐based management. The need to move beyond models of single stressor effects on ecosystems and a reliance on setting universal environmental limits has been difficult to operationalize and translate into management actions. This is especially important in multi‐use estuarine and coastal ecosystems where stressors are not spatially separated. New opportunities to develop a mechanistic understanding of how stressors interact will allow us to develop and test frameworks for addressing cumulative effects that go beyond additive models of multiple stressors, which despite wide recognition of the problem still dominate the literature (Hodgson and Halpern 2019, O'Brien et al. 2019). Cumulative risk assessment frameworks are being trialed but at the heart of these models is the need to understand the mechanism linking stressors to ecosystem consequences (Hodgson and Halpern 2019, Stelzenmüller et al. 2020). Managers and policy makers are often charged by regional, national policy and international obligation to try to navigate towards more positive futures by reversing the trajectories of biodiversity loss, diminution of ecosystem services for coastal and estuarine ecosystems, but operationalizing actions to tackle multiple stressors remains a frontier. Knowledge of interaction networks between multiple environmental drivers and ecosystem components extends the prospect of positive interventions not only through managing stressor loads but also by enhancing or restoring ecosystem responses that mitigate adverse effects.

Empirical validation of the potential for regime shifts and tipping points in ecological systems justifies to managers and society at large that these phenomena are real, common and serious. The effects we observed in EINs associated with sediment and nutrient pollution highlight the need for multiple complementary management actions and a move away from approaches that focus on stressors and not ecosystem responses. Insular actions on individual stressors based on assumptions of gradual change, involving the application of “set and forget” environmental policies risk unintended consequences (Thrush et al. 2016). More integrative ecosystem‐based management requires empirical research working alongside theoretical ecology to create new opportunities for managers to develop new forms of risk assessment and actions to manage for surprise.

Appendix S1

Click here for additional data file.