Rising tides, cumulative impacts and cascading changes to estuarine ecosystem functions

In coastal ecosystems, climate change affects multiple environmental factors, yet most predictive models are based on simple cause-and-effect relationships. Multiple stressor scenarios are difficult to predict because they can create a ripple effect through networked ecosystem functions. Estuarine ecosystem function relies on an interconnected network of physical and biological processes. Estuarine habitats play critical roles in service provision and represent global hotspots for organic matter processing, nutrient cycling and primary production. Within these systems, we predicted functional changes in the impacts of land-based stressors, mediated by changing light climate and sediment permeability. Our in-situ field experiment manipulated sea level, nutrient supply, and mud content. We used these stressors to determine how interacting environmental stressors influence ecosystem function and compared results with data collected along elevation gradients to substitute space for time. We show non-linear, multi-stressor effects deconstruct networks governing ecosystem function. Sea level rise altered nutrient processing and impacted broader estuarine services ameliorating nutrient and sediment pollution. Our experiment demonstrates how the relationships between nutrient processing and biological/physical controls degrade with environmental stress. Our results emphasise the importance of moving beyond simple physically-forced relationships to assess consequences of climate change in the context of ecosystem interactions and multiple stressors.

As we continue to alter coastal ecosystems, the ability of estuaries to deliver multiple ecosystem services decreases, but our need for the benefits they confer grows. This discrepancy between our reliance on estuaries and their capacity to mitigate change can lead to ecosystem collapse [1][2][3][4] . Estuaries play vital roles in regulating nutrient flux. Some estuaries have been shown to reduce nitrogen (N) stocks by over 70% through denitrification alone, mitigating anthropogenic alterations of coastal ecosystems [5][6][7] . However, estuaries bear the brunt of human induced stress. As the transition between land and sea, estuaries are subject to both terrestrial and marine stressors that can work synergistically to reduce overall functionality and productivity. Climate change (including sea level rise, increased wave action, changes in precipitation patterns, warming, etc.) and intensification of land use promote sediment and nutrient runoff [8][9][10][11][12] . The combination of stressors acting simultaneously constitutes a multiple stressor scenario where predicting the consequences of environmental change on ecosystem function requires examination of the interactions between physical and biological processes that create ecosystem networks. More simplistic single stressor cause-and-effect relationships are likely of limited value because estuarine function is based on complex interactions. Although complexity poses major challenges, we must exploit knowledge of key connections between processes to develop comprehensive models of interaction networks.
Microphytobenthos (MPB) are often the dominant source of primary production on estuarine intertidal flats and play important roles in ecosystem function [13][14][15][16][17] . Shifts from benthic to planktonic production are associated with eutrophication and/or increased turbidity. In estuaries where water column nutrients limit primary production, the MPB regulate ammonia flux across the sediment water interface [17][18][19] . With their rapid growth, high turnover rate, and high palatability, MPB respond quickly to nutrient pulses and aid in the transfer of nutrients to higher trophic levels for storage, transport and further processing. MPB have other important functional roles in estuarine habitats. The coagulation of surface sediments by mucus produced by MPB reduces nutrient release from buried substrate and regulates turbidity by inhibiting the resuspension of particulates 18 . In healthy estuaries, these regulatory processes are tightly coupled 14,20 . Healthy ecosystems can cycle excess nutrients through microbial processing, such as denitrification [21][22][23] . However, in disturbed systems, ecosystem processes may be decoupled. Measuring any single parameter may not provide a complete picture. Therefore, understanding the interactions between ecosystem processes and biological or physical controls provides a more comprehensive descriptor of ecosystem health.
Informed by knowledge of ecosystem processes, we experimentally assessed changes in biological and physical controls of nutrient flux. This involved the deployment of in-situ mesocosms (area = 1 m 2 , volume = 180 L) in an experimental design that manipulated sea level (+18 cm), nutrient content (+87 g N m −2 , 7 g P m −2 ), and mud content (5 mm deposition event) individually and in combination to mimic realistic anthropogenic disturbance on intertidal flats. Experimental results were compared with processes measured along elevation gradients as a method of substituting space for time and assessing the validity of our manipulations. We measured nutrient flux, sediment characteristics, and microphytobenthic production (estimated by benthic chlorophyll-α content) to quantify changes in ecosystem function. We expected a strong correlation between site characteristics and nutrient flux in unaltered (control) experimental plots and transect sites. In contrast, sites altered by increased sea level (SLR), nutrient content, and/or sediment content, were expected to display shifts in the relationships between site characteristics and flux.

Results
Benthic production. Benthic production was increased by nutrient additions, decreased by sediment additions, and unchanged when both sediment and nutrients were added (Fig. 1). Sea level had the greatest effect on benthic chlorophyll-α content, and masked any effect of nutrient and sediment additions (Fig. 1). Transect results were consistent with experimental plot data and indicated a decline in MPB production with decreasing elevation/increased inundation (R 2 = 0.79, p < 0.01, n = 54).
Modelling ecosystem networks. Generalized linear models were developed using both experimental plot and transect data to determine which interactions were most important in regulating nutrient flux (our proxy for ecosystem function). Models developed of ammonia flux from experimental plots were consistent with transect data and indicated that factors controlling ammonia flux differed between incubations conducted in the light (photosynthetically active, equation 1) and incubations conducted in the dark (photosynthetically inactive, equation 2). Models were developed using all the environmental variables assessed in this study: nutrient fluxes, organic matter content, porosity, O 2 consumption, and benthic production. However, the most parsimonious models were: . where PO 4 3− is the flux of soluble reactive phosphorus in ug m −2 hr −1 , chla is the concentration of benthic chlorophyll-α in mg m −2 , O 2 is oxygen consumption in mg m −2 hr −1 , PR is photosynthetic rate in g-Carbon (mg chlorophyll-α * hr) −1 , ρ is porosity (fraction), and NO x − is the flux of combined nitrate and nitrite in ug N m −2 hr −1 . The flux of ammonium (NH 4 + ), under both light and dark conditions, was tightly coupled to soluble reactive phosphorus (SRP) flux and decreased with decreasing elevation. Ammonia flux in light conditions was dominated by primary productivity with chlorophyll-α, O 2 consumption, and photosynthetic rate as regulatory parameters. Ammonia flux from dark chambers was dominated by respiration, including porosity, which affects organism movement and diffusive flux of NO x − produced through microbial processing, and O 2 consumption. The models for both light and dark NH 4 + were most successful at describing plots with no stressors, and model fit decreased as the number of stressors increased (Fig. 2). The decline in R 2 occurred regardless of stressor or combination of stressors within each plot.

Discussion
Our experimental plots showed a positive correlation between MPB production and NH 4 + flux, demonstrating the interacting links between MPB, macrofauna, and nutrient flux. MPB regulate NH 4 + flux by capping sediments with mucus, direct uptake, and the provision of food resources to higher trophic levels [24][25][26] . The tight coupling between SRP and NH 4 + fluxes retained across all experimental treatments is indicative of remineralisation and release across the sediment water interface (Fig. 3a). The observed relationship between porosity in dark chambers and benthic production in light chambers indicates that MPB do play a role in regulating the exchange of nutrients across the sediment water interface in our experimental plots. Benthic grazers also play a role in nutrient flux as they secrete NH 4 + as waste and bioturbate, releasing stored nutrients and oxygenating sediments (Fig. 3a) 27 . Increased food stores attract macrofauna, further increasing bioturbation and the subsequent release of nutrients and greater O 2 consumption [28][29][30][31] . Therefore, an increase in MPB stocks supports benthic grazers and leads to a rise in the release of NH 4 + and SRP from sediments 27,31,32 . In our experimental plots, we observed a tight coupling between SRP and O 2 consumption which is an indicator of bioturbation and bioirrigation. Bioturbating macrofauna consume O 2 through respiration and their disruption of the anoxic microzone stimulates microbial processing and further O 2 consumption. The increase in NH 4 + flux observed with increased MPB biomass and greater photosynthetic rate emphasise both the importance of sediment nutrient cycling to primary production in non-eutrophic systems, and the roles of bioturbation and bioirrigation in elevating solute transport rates.
As biogeochemical gradients are altered by benthic grazers, microbial processing of nutrients and organic matter are stimulated. In particular, anammox, coupled nitrification/denitrification, and respiration process nutrients and organic matter released from bioturbation [33][34][35][36] . In dark conditions, NO x − flux replaced MPB biomass and photosynthetic rate as regulating factors of NH 4 + flux. NH 4 + and NO x − are directly released into the  water column if the nutrients are not utilized at the sediment water interface. Moreover, nutrients are more readily released from porous sediments, further boosting the flux of both NH 4 + and NO x − , evidenced by the importance of porosity for dark nutrient flux. In addition to direct release, NO x − can be produced through nitrification, which converts NH 4 + to NO 3 − , using O 2 as an electron acceptor. Because benthic O 2 consumption increases as NO x − flux increases in our experiment, this may indicate nitrification occurring within the system. The net result of these changes in flux indicate these systems can shift from a sink of bioavailable nutrients to a source as the interaction network degrades. The models we have developed as a result of our experimental data are informed by knowledge of ecosystem process and illustrate how subtle changes in environmental conditions can shift ecosystems in profound ways (Fig. 3).
The loss of MPB production with SLR in experimental plots and decreased elevation along transects may be linked to a rise in planktonic primary production. In ecosystems reliant on MPB production, a shift from benthic to planktonic production can result in a positive feedback loop. As nutrients are added to the system, primary production shifts to the water column (Fig. 3b) 37,38 . With the transfer of production from the benthos to the water column, the sediments are shaded and benthic primary production decreases, resulting in increased nutrient release, sediment mobility, turbidity, and eutrophication ( Fig. 3b) 1,13,39,40 . Environmental change such as sediment loading and SLR, can thus contribute to the movement of primary production from benthic to planktonic. Suspended sediments increase turbidity and smother benthic algae 8,11,25,27 . Similarly, when intertidal habitat is converted to subtidal, the amount of light available decreases, potentially shifting a higher proportion of production to the water column. This results in the loss of MPB control over nutrient flux across the sediment water interface and can lead to a positive feedback loop fueling eutrophication and a significant shift in the interaction network (Fig. 3b) 41 .
Results show that our ability to predict NH 4 + flux is compromised as the number of stressors increases, which indicates a breakdown in the networked processes governing ecosystem function. Our experiments also indicated that nutrient and sediment additions effects, although very subtle at the low levels used in our experiment, were masked by SLR. The non-linear, non-additive effects of nutrient pollution, sediment addition, and SLR indicate that model predictions of future change will need to weight certain stressors more heavily. Currently, there is no clear way of predicting a tipping point 3,4 . However, one useful approach is to assess the risk of change to environmental stress in terms of loss of feedbacks and networked interactions 10,14,[42][43][44] . Based on our experiment and survey of patterns apparent along gently sloping intertidal to subtidal sandflats, we conclude that SLR is the most important factor in determining an environmental tipping point for primary producers. In this experiment, an 18 cm increase in SLR decreased benthic chlorophyll-α concentrations by ~30%. With SLR estimates ranging from 18-59 cm in New Zealand by 2090 45 , losses in MPB biomass could be higher. As we continue to alter coastal ecosystems, we will observe varied and often unexpected environmental responses. Not only will individual parameters change, but interactions or relationships between these factors will also be altered. We cannot rely on single descriptors to identify tipping points. Instead, we should characterise ecosystem networks to identify at risk communities and potential solutions.
At each site, 3 replicates of each of the following treatments were established: control (no amendments/additions), +nutrients, +sediment, simulated SLR, +nutrients and sediment, +nutrients and SLR, +sediment and SLR, and all three stressors combined (Fig. S1). Once treatments were established, plots settled for 5 months before sampling occurred. Benthic chlorophyll-α content, porosity, and sediment organic matter content were collected during deployment and again in March 2016. Nutrient fluxes were measured in spring (November) 2015 and late summer (March) 2016 to assess flux of NO x − , SRP, NH 4 + , and O 2 across the sediment water interface. To compare results of in-situ SLR manipulations, transects (length = 100 m, sampled every 20 m) were established perpendicular to shore just above the study site and extending to the subtidal zone (Fig. S1). This allowed us to compensate for any artefacts of in-situ manipulation of SLR.

Mimicking environmental change.
Estimates of SLR for New Zealand range from 18-59 cm by 2090 45 .
To be conservative, we use a simulated SLR of 18 cm. Cylinders (diameter = 1.60 m, height = 0.36 m) were constructed using opaque white HDPE plastic sheets. Cylinders were pushed into sediment until 18 cm of plastic remained above the sediment surface (Fig. S2) to create mesocosms. Due to the sediment grain size and texture, the mesocosm sealed against the sediment, retaining water. Water was exchanged at high tide, but remained subtidal at low tide (Fig. S2). Cylinders were installed to represent 4 different treatments: SLR, SLR + sediments, SLR + nutrients, SLR + sediment + nutrients (Fig. S1 and S2). For treatments with added nutrients, 450 g of controlled release fertiliser (Osmocote Total All Purpose) was inserted 15 cm below the sediment surface using a sediment corer (87 g N m −2 , 7 g P m −2 ). Nutrient dosage and methodology was adapted from Douglas et al. 2016 47 who performed experiments in similar systems. Sediment additions (adapted from Lohrer et al. 2004) 8 were achieved using slurries of marine mud in seawater. Sediments were collected from the field, sieved to remove shell hash and other large particles, and allowed to sit in fresh water for 2 weeks to neutralise any macrofauna. Sediments were added at low tide and spread evenly across plots to simulate a 5 mm deposition event.
Scientific RepORTS | 7: 10218 | DOI:10.1038/s41598-017-11058-7 Sediment Characterization. Sediment samples were collected using a 20 mL syringe to 2.5 cm depth to keep a constant wet volume and dried at 60 °C for 2 days. Porosity was determined using equation (3): where P is porosity, V is volume of wet sediment, W is weight of wet sediment, and D is weight of dry sediment. SOM percentages were measured using loss on ignition 48 . Samples were dried at 60 °C for 2 days and combusted at 525 °C for 4 hours. SOM content was determined using equation (4): Nutrient flux. Concentrations of NH 4 + , NO x − , and SRP were collected twice (November 2015, March 2016) during a midday high tide to maximize photosynthetic capacity. Chambers (diameter = 15 cm, volume = 1 L) were placed over the sediment surface to determine fluxes over approximately 6 hours (Fig. S3). Light and dark chambers were deployed side by side within each plot. Light chambers were constructed of translucent plastic to allow photosynthetic activity to continue. Dark chambers were constructed of black plastic painted white to reduce thermal gradients but stop photosynthetic activity. Flux rates were determined using equation (5).
F I where C I and C F are the initial and final concentrations of an analyte, respectively, T is incubation time (h), and A is the surface area of the core (m 2 ) 49 . Initial samples were collected just after deployment and final samples were collected after approximately 6 hours. To account for water column production, ambient samples were collected along with initial samples and water was incubated in light and dark bottles simultaneously with chamber samples. Water samples were filtered using Whatman GF/F filters (pore size of 0.7 µm) and analysed with a Lachat Quick-Chem 8000 automated ion analyser for NO x

Microphytobenthic production.
Chlorophyll-α concentrations were used as a proxy for MPB production. Samples were collected just after deployment and again in March 2016. Sediment cores (area = 1.13 cm 2 , depth = 1 cm) were collected in triplicate within each experimental plot and transect site (6 sites per transect). Samples were frozen immediately and processed within 1 month of collection. Chlorophyll-α was extracted from sediments for approximately 18 hours at 0 °C in a solvent mixture of 45:45:10% methanol: acetone: deionised water 21 . After extraction, samples were vigorously mixed and allowed to settle before analysis using a Shimadzu spectrophotometer 50,51 . Samples were acidified to account for phaeophytin concentrations.

Statistics.
A generalized linear model (GLM) was developed to investigate the individual impacts of each stressor, the combined effects, and determine the relative importance of each. Linear regressions were used to assess trends in benthic production with elevation for transect sites. To determine differences in chlorophyll-α content between groups, a Dunn's Test was conducted. All data were analyzed using R. Errors reported are standard errors.