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.

This dataset contains the data used for the Ecological Applications publication "Cumulative stressors reduce the self-regulating capacity of coastal ecosystems".  It contains two sheets, labelled as per the use in the manuscript.

Study design

Using local knowledge, we identified 24 sites in 15 estuaries across New Zealand (Figure 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 (see methods below) 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).

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 m² plots for each of 3 treatments, high fertilizer 600 g N m^-2, medium fertilizer 150 g N m^-2, and disturbance controls.  We used Nutricote® 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 min (fixed vertically 10 cm above sediment surface; Odyssey, Dataflow Systems Pty Ltd, Christchurch, NZ) and sediment temperature sensors logging temperature every 30 min (buried at depths of 3 and 7 cm; i-Buttons, Maxim Integrated Products Inc., 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 dia, 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 (particle size range 0.01 to 3500 µ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 10-AU fluorimeter (Arar and Collins 1997).  Sediment organic matter content was determined by weight loss on ignition (60^oC for 24 h or until stable mass + combustion for 4 h at 550^oC) (Dean 1974).

Two macrofauna cores (13 cm dia, 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 greater than 20mm 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) were deployed in chambers in case temperature was required for the EINs. Biochemical measures made were net sediment-water O₂, N₂, and NH₄ fluxes. NO_x 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 cm x 50 cm x 15 cm height) pressed 5 cm into the sediment during low tide.  On the incoming tide (water depth ~50 cm), Perspex domes sealed ~40 liters 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 x 60 ml syringe for solute, and 2 x 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 O₂ concentrations in the water samples were measured using an optical probe (ProODO YSI Inc.). Samples were then filtered through a 0.45 µm Whatman GF/C glass fiber filter and stored frozen (-20 °C) prior to analysis of NH₄⁺-N Water from airtight syringes were preserved in zinc chloride and stored in gas tight exetainers (Labco, UK) at 4°C until analysis of N₂ concentration. Solute and gas fluxes were calculated as (C_end-C_initial * V ) / A * T, where C is nutrient or oxygen concentration (µM L^-1), V is the volume of seawater inside the chamber (L), A is the area of sediment enclosed by the chamber (m²), and T is the elapsed time between initial and final samplings (h). 

Water samples collected from the chambers were analyzed for NH₄⁺ with a Lachat QuickChem 8000 Series FIA+ (Zellweger Analytics Inc. Milwaukee, Wisconsin, 53218, USA) using standard operating procedures for flow injection analysis (NH₄⁺-N detection limit of 0.07 mmol L^-1).  N₂ gas fluxes were analyzed on a quadrupole membrane inlet mass spectrometer (MIMS; Bay Instruments, Cambridge MD, USA) using the N₂/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 non-linear 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, NH₄ 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 non-linear response). Linear regressions detected no significant relationships and explained less than 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 non-linear and the breakpoints detected were useful.

Regression tree analyses were conducted using the RPART package (Therneau T. 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²) by ≥0.03 and each terminal node (final most homogeneous group) had to contain at least 10 data points (i.e., 2 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 450 μM m^-2 s^-1 PAR (Figure 3). Once we knew that the regression tree analysis not only suggested a non-linear 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 was broken into three groups: Clear sites with >450 μM m^-2 S^-1 PAR, turbid sites with <350 μM m^-2 S^-1 PAR; and an intermediate set with too few data points to warrant further analysis.

To use this data, you need to understand marine sediment ecology and the biogeochemisty assocoated with nutrient and sediment fluxes.