DOI: 10.5061/dryad.8cz8w9gms
title: Cumulative stressors reduce the self-regulating capacity of coastal ecosystems

description: 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.
description: <strong>Study design</strong>
 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 &lt;0.25). Sites were
 situated at about  mid-tide level, on sandy permeable sediments and in
 areas occupied by the functionally important bivalves <i>Macomona
 liliana</i> and <i>Austrovenus stutchburyi</i> at their
 natural densities (see Results). At each site we established three 9
 m<sup>2</sup> plots for each of 3 treatments, high fertilizer
 600 g N m<sup>-2</sup>, medium fertilizer 150 g N
 m<sup>-2</sup>, 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. <strong>Data
 collection</strong> 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<sup>o</sup>C for 24 h or until stable
 mass + combustion for 4 h at 550<sup>o</sup>C) (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. <i>Macomona</i> and <i>Austrovenus</i>
 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<sub>2</sub>, N<sub>2</sub>, and
 NH<sub>4</sub> fluxes. NO<sub>x</sub> 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<sub>2</sub> 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<sub>4</sub><sup>+</sup>-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<sub>2</sub>
 concentration. Solute and gas fluxes were calculated as
 (<i>C<sub>end</sub>-C<sub>initial</sub> *
 V</i> ) / <i>A</i> * <i>T</i>, where
 <i>C</i> is nutrient or oxygen concentration (µM
 L<sup>-1</sup>), <i>V</i> is the volume of
 seawater inside the chamber (L), <i>A</i> is the area of
 sediment enclosed by the chamber (m<sup>2</sup>), and
 <i>T</i> is the elapsed time between initial and final
 samplings (h).  Water samples collected from the
 chambers were analyzed for NH<sub>4</sub><sup>+
 </sup>with a Lachat QuickChem 8000 Series FIA+ (Zellweger Analytics
 Inc. Milwaukee, Wisconsin, 53218, USA) using standard operating procedures
 for flow injection analysis
 (NH<sub>4</sub><sup>+</sup>-N detection limit of
 0.07 mmol L<sup>-1</sup>).  N<sub>2</sub> gas
 fluxes were analyzed on a quadrupole membrane inlet mass spectrometer
 (MIMS; Bay Instruments, Cambridge MD, USA) using the
 N<sub>2</sub>/Ar method (analytical precision &lt;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.
 <strong>EIN analyses</strong>
 While we had <i>a priori</i> 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<sub>4</sub> 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<sup>2</sup>) 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<sup>-2 </sup>s<sup>-1</sup> 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 &gt;450 μM m<sup>-2
 </sup>S<sup>-1 </sup>PAR, turbid sites with &lt;350
 μM m<sup>-2 </sup>S<sup>-1</sup> PAR; and an
 intermediate set with too few data points to warrant further
 analysis.
description: To use this data, you need to understand marine sediment ecology
 and the biogeochemisty assocoated with nutrient and sediment
 fluxes.