The Effects of Long-Term Nitrogen Enrichment on Estuarine Benthic-Pelagic Coupling

: Biogeochemical cycling in the marine coastal zone regulates the availability of nitrogen and carbon within soft sediment habitats. However, these pathways are being fundamentally altered by anthropogenic increases in nutrient delivery. Few studies have incorporated long-term enrichment and ecological complexity (in situ experiments), restricting our ability to manage effectively and prevent ecological shifts. This study investigates the inﬂuence of sediment nutrient availability (at 3 levels, across 2 seasons) on biogeochemical cycling over a 20-month period in 4 estuaries. Overall, net denitrification rates were highly variable, ranging between 4 and 208 µ mol N m − 2 h − 1 . However, no increases were observed with increasing enrichment highlighting the limited capacity for nitrogen removal in response to large increases in bioavailable nitrogen. Additionally, macrofaunal communities and sediment trophic status were shown to have important inﬂuences on nitrogen processing. Overall, alterations to ecosystem relationships and the appearance of non-linear responses to increasing nutrient enrichment reveal the vulnerability of estuaries to increasing stressor loads owing to the increased likelihood of reaching a tipping point.


Introduction
Benthic-pelagic coupling within soft sediment coastal ecosystems regulates the transfer of mass and energy via biogeochemical transformations and biological interactions, which ultimately determines the production, structure and stability of both systems [1][2][3].An important component of benthic-pelagic coupling is nutrient cycling, which alters the uptake, transformation and removal of bioavailable nitrogen (N) as well as regulating the production and metabolism of carbon (C) [4,5].The processes and pathways involved in both N and C cycling are increasingly affected by a multitude of pressures facing coastal ecosystems.In particular, the anthropic acceleration of the global N cycle has resulted in sources exceeding sinks by more than 40% [6], due to an exponential increase in agricultural productivity leading to the export of N to coastal ecosystems [7,8].If this N enrichment exceeds the capacity for assimilated-enhanced benthic production and/or removal then a system can be pushed beyond a tipping point whereby pelagic primary production dominates and coastal eutrophication persists [9,10].
Responses of oligotrophic unvegetated, soft sediment habitats to increasing N enrichment include the direct stimulation of benthic microalgae [11] which can result in increases in the abundance and biomass of macrofaunal grazing populations [12].Changes to benthic microalgae production, macrofaunal diversity and their interactions alters oxygen distribution, sediment biogeochemistry and organic matter (OM) availability and therefore can directly affect sediment nutrient processing [11,13].Sediment nutrient processing is a globally valuable ecosystem function, influencing the supply and flux of nutrients both within and between marine habitats [14].Of particular importance is sediment denitrification (DNF), a process that can remove a significant proportion of excess bioavailable N (up to 80%) [15] and therefore has been suggested to aid in the resilience of coastal ecosystems to eutrophication [16].
DNF is a heterotrophic, anaerobic process whereby nitrate is used as an electron acceptor to oxidise OM [17].Owing to its considerable importance, the factors controlling DNF are reasonably well understood and are largely attributed to the availability of oxygen, nitrate and carbon (through the decomposition of OM) [5].OM, nitrate and oxygen are in turn regulated by ambient concentrations in the overlying water column (dependent on local hydrodynamic conditions) and within the sediment by benthic primary producers, macrofaunal activities and sediment properties [17][18][19][20].The positive influence of macrofaunal activities on DNF and nutrient processing is however, likely to be diminished if communities are negatively affected by increasing nutrient enrichment [13].This can occur through ammonium toxicity, the formation of hydrogen sulphide and/or increased anoxia within the sediment [21][22][23].While the response of macrofaunal communities to sediment N enrichment is species and location dependent [23,24], understanding the response of key macrofaunal species that disproportionally contribute to ecosystem function is likely to be a crucial component in understanding any alterations in coastal biogeochemical cycling.
The delivery of N to the coastal zone is intensifying, invoking changes to ecosystem interactions and increasing the likelihood of nonlinear responses, which can collectively degrade ecosystem functionality and reduce the ability of an ecosystem to adapt to further stress [25,26].It is therefore useful to incorporate real world ecological complexity through in situ investigations, in order to predict broadscale ecosystem responses [27].For example, in situ nutrient enrichment experiments have increased our understanding of the biological and biogeochemical changes that occur with eutrophication in unvegetated sediments [13,27,28].However, most are restricted both spatially (to a few locations) and temporally (to a duration of a few weeks or months [29][30][31].This negates the incorporation of longer-term process such as macrofauna recruitment cycles and shifts in community structure as the system equilibrates to the added stressor, decreasing the real-world relevance.Moreover, the spatially restricted nature of many experimental studies hinders quantifying the natural variability in the response, reducing our ability to generalise and inform management. This study builds upon a preceding in situ experiment by Thrush, et al. [27] where sediment nitrogen concentration was enriched across 15 estuaries over a 9-month period, demonstrating that the interaction of multiple stressors (nutrient enrichment and water column turbidity) can cause distinct changes in ecosystem interaction networks that drive nutrient processing.This study extended the duration of nutrient enrichment (to 20 months) in a restricted number of estuaries (4) to investigate how long-term in situ nutrient enrichment alters nutrient processing and biogeochemical cycling.By integrating over a longer enrichment period and including additional measures of biogeochemical cycling recorded in two seasons, the present study aims to quantify the variability in response, which is essential for prediction and therefore an integral component to aiding the management of multiple stressors within coastal ecosystems.We hypothesise that high nutrient enrichment will negatively affect the abundance of key macrofaunal species, which in turn will influence nutrient processing.Additionally, we hypothesise that DNF rates will increase with nutrient enrichment before plateauing at high enrichment.

Study Sites
This study continued a large scale field experiment described in Thrush, et al. [27], where sediments were enriched with nitrogen for 9 months at 24 sites (in 15 estuaries) which encompassed a gradient in water column turbidity.Whereas Thrush, et al. [27] investigated alterations in ecological interaction networks, this study selected a reduced number of sites (6 sites in 4 estuaries) to identify treatment effects on nutrient processing over a longer time period (20 months).This extended period ensured the incorporation of longer-term processes such as any alterations to macrofaunal communities and sediment properties.We also included measurements in winter as well as early summer in order to better resolve seasonal impacts.
The 6 sites were located within the North Island of New Zealand within 4 shallow, barrier enclosed estuaries (Figure 1).Each estuary had an extensive intertidal area and semi-diurnal tides.Sites encompassed natural spatial heterogeneity in the abundance of two functionally important bivalve species (the suspension-feeding clam Austrovenus stutchburyi and deposit-feeding wedge shell Macomona liliana) and sediment properties (organic and mud content), both of which can substantially alter biogeochemical processes [32,33].All of the sites were positioned within the mid-intertidal region and located along a gradient in water column turbidity such that light availability (measured as photosynthetically active radiation (PAR)) reaching the sediment surface ranged from a median of 38-283 µmol m −2 s −1 during tidal submergence [34].

Study Sites
This study continued a large scale field experiment described in Thrush, et al. [27], where sediments were enriched with nitrogen for 9 months at 24 sites (in 15 estuaries) which encompassed a gradient in water column turbidity.Whereas Thrush, et al. [27] investigated alterations in ecological interaction networks, this study selected a reduced number of sites (6 sites in 4 estuaries) to identify treatment effects on nutrient processing over a longer time period (20 months).This extended period ensured the incorporation of longer-term processes such as any alterations to macrofaunal communities and sediment properties.We also included measurements in winter as well as early summer in order to better resolve seasonal impacts.
The 6 sites were located within the North Island of New Zealand within 4 shallow, barrier enclosed estuaries (Figure 1).Each estuary had an extensive intertidal area and semi-diurnal tides.Sites encompassed natural spatial heterogeneity in the abundance of two functionally important bivalve species (the suspension-feeding clam Austrovenus stutchburyi and deposit-feeding wedge shell Macomona liliana) and sediment properties (organic and mud content), both of which can substantially alter biogeochemical processes [32,33].All of the sites were positioned within the mid-intertidal region and located along a gradient in water column turbidity such that light availability (measured as photosynthetically active radiation (PAR)) reaching the sediment surface ranged from a median of 38-283 µmol m −2 s −1 during tidal submergence [34].

Experimental Design
The allocation of nitrogen, plot size and level of replication within each site was determined by the experiment of Thrush, et al. [27].At each site, sediment nutrient concentrations were elevated at two levels resulting in three treatments (control: 0; medium: 150; high: 600 g N m −2 ) across nine (n = 3) 9 m 2 experimental plots.The uniform elevation of porewater N in each plot was achieved following Douglas, et al. [35], where a sediment plug (3 cm diameter, 15 cm deep) was extracted from 180 evenly spaced holes using a handheld corer, a known volume of slow-release urea fertiliser (Nutricote; 40:0:0 N:P:K) was added and then the sediment plug replaced.Urea fertiliser was chosen as it quickly hydrolyses to ammonium (NH 4 + ), a product released during the remineralisation of OM and a typically limiting nutrient within New Zealand estuaries [36].In addition, these levels of N enrichment have been shown to alter benthic macrofaunal diversity but do not result in defaunation [13].While this method of N enrichment does not necessarily mimic the process of increased OM degradation, it does allow the stimulation of porewater enrichment associated with long-term eutrophication [32,35].To maintain nutrient enrichment over the full experimental period, fertiliser was applied in March and November 2017 and following the sampling in June 2018.This resulted in a minimum of 5 months between the application of fertiliser and sampling (Figure 1).
Sampling first took place after 15 months of sediment nutrient enrichment in June 2018 (winter; 3 sites TAU-T, TAU-O, RAG) and after 20 months in November 2018 (early summer; all 6 sites), with care taken not to resample the same area within each plot.Data collected by Thrush, et al. [27] in November 2017 was not included in this study as the full suite of measurements made here were not available.The period of sediment enrichment before sampling provided the opportunity to incorporate longer-term processes, while the seasonal component examined temporal variations.

Field Sampling
At three sites during winter (TAU-O, TAU-T and RAG) and all sites during summer, benthic incubation chambers were deployed over a full submerged tidal period of approximately 4 h to assess solute and dissolved gas fluxes.During the preceding low tide, two chamber bases (50 cm × 50 cm × 15 cm) per plot were inserted 5 cm into the sediment and equipped with a light and temperature logger (HOBO Pendant ® , Onset Computer Corporation, Bourne, MA, USA) and a dissolved oxygen logger (PME miniDOT).On the incoming tide, a Perspex dome lid (enclosing ~40 L) was placed onto each chamber base and all air bubbles were removed before each chamber was sealed.An opaque shade cloth covered one chamber lid per pair (within each plot) to prevent exposure to sunlight, while the other was left to ambient light conditions.Once sealed and at the end of the incubation, seawater samples (1 × 60 mL syringe for solute, 2 × 60 mL airtight syringes for gas concentrations) were extracted from each chamber.Solute samples were filtered on-site through a 0.45 µm Whatman GF/C glass fibre filter and frozen at −20 • C until analysis, and gas samples were transferred to airtight exetainers (Labco, UK) and preserved with zinc chloride at 4 • C until analysis.
To assess the vertical variation in sediment OM degradation within each plot, rapid organic matter assay (ROMA) plates were deployed between 7 and 10 d prior to the benthic chamber measurements during the summer sampling, following the methods described by O'Meara, et al. [37].In brief, each plate consisted of 3 × 0.9 mL wells at depths of 1, 3, 5, 7, 10 and 15 cm, and filled with an agar mixture consisting of 0.029 g C mL −1 .Once removed, each plate was gently washed to remove any sediment before measuring the change in agar volume per well within 6 h of collection.
For determination of sediment characteristics (sediment chlorophyll a, phaeopigments, grain size and organic content) during both sampling events, five pooled sediment cores (2.6 cm dia) were collected from each experimental plot at a depth of 0-2 cm and frozen at −20 • C until analysis.Three separate replicates of four pooled sediment cores were additionally taken at depths of 0-2 cm and 5-7 cm for analysis of sediment porosity, porewater dissolved inorganic nutrient concentrations and fluorescence characteristics of dissolved organic matter (DOM).One macrofauna core (13 cm dia, 15 cm depth) was taken from each plot (3 per treatment per site) and sieved on a 500 µm mesh before being preserved in 70% Isopropyl alcohol.Holes left by the macrofauna cores were filled with clean sand to prevent erosion/collapse of the surrounding sediment.

Laboratory Analysis
To determine sediment properties, samples were thawed and homogenised.Organic content was determined by weight loss on ignition, where samples were dried at 60 • C until a constant weight and then combusted at 550 • C for 4 h.For grain size determination, samples were digested in 10% hydrogen peroxide before being measured by laser diffraction (Malvern Mastersizer-3000).Chlorophyll a and phaeopigments were extracted from freezedried sediment using 90% buffered acetone and then measured before and after acidification using a fluorometer (Turner Designs 10-AU).
Samples for sediment porosity and porewater extraction were processed within 24 h of collection.To calculate sediment porosity, the difference in wet and dry weight was determined after drying (60 • C) for 7 d or until a constant weight.To extract sediment porewater for determination of dissolved inorganic nutrient concentrations and fluorescence characteristics of DOM, 4 mL of de-ionised water was first added to each sample, before being vortexed, centrifuged and then filtered through a 0.45 µm Whatman GF/C glass fibre filter and stored at −20 • C. Analysis of dissolved inorganic nutrient concentrations (from the chamber incubations and extracted porewater; NH 4 + -N, NO 3 − -N, NO 2 − -N and PO 4 3− -P) was conducted using standard operating procedures for flow injection analysis on a Lachat QuickChem 8000 Series FIA+ (Zellweger Analytics Inc.).Fluorescence characteristics of porewater DOM were analysed using a fluorescence spectrometer (Jobin Yvon Aqualog ® , Horiba, Kyoto, Japan) with a 2 s integration time, 3 nm step-size and a measurement range of 240-600 nm excitation and 245-800 nm emission.
Analysis of gas samples taken from the dark chambers were performed using Membrane Inlet Mass Spectrometry (MIMS) as per O'Meara, et al. [38].Dark only measurements were chosen to exclude the influence of photosynthetic O 2 production by microphytobenthos (MPB), the dominant primary producer at each site [34], as increases in oxygen saturation can disrupt N 2 quantification for assessment of DNF rate [39].Transformation of raw data into dissolved gas concentration (µmol L −1 ) was based on gas solubility, which was calculated using in situ temperature, pressure and salinity measurements (taken at the same time of sampling) as per Hamme and Emerson [40].Variations in N 2 concentrations were normalised with Argon (Ar), a biologically inert gas, and the resulting N 2 /Ar ratios used to calculate N 2 gas concentration.
Macrofauna samples from the summer sampling (3 per treatment per site) were stained with Rose Bengal and fauna separated from any shell hash before being identified to the lowest possible taxonomic level (usually species) and counted.For the winter samples, only adult (>10 mm) Austrovenus stutchburyi and Macomona liliana were counted.

Parameter Derivations
Fluxes of dissolved oxygen, inorganic nutrients and N 2 gas were calculated as the difference between final and initial concentrations, before being standardised by chamber volume, sediment surface area and incubation time.Oxygen fluxes from the dark chambers are presented as sediment oxygen consumption (SOC; µmol O 2 m −2 h −1 ), where a positive net flux represents uptake by the sediment.Oxygen fluxes from the light chambers represent net primary production (NPP; µmol O 2 m −2 h −1 ) and are not presented individually but have been previously described by Mangan, et al. [41] and are used to calculate trophic status (see below).A net positive or negative nutrient flux represents an efflux from and influx into the sediment, respectively.Fluxes of NO 2 − , NO 3 − and PO 4 3− were close to detection limits and therefore were not presented individually but NO 2 − and NO 3 − were included in the calculation of total dissolved inorganic nitrogen flux.N 2 fluxes represent a balance between N fixation and production, hereafter DNF, and therefore a positive net flux was attributed to net DNF, produced via DNF and/or anammox pathways.Negative fluxes (3/81) were omitted based on the introduction of bubbles within the sample.
In the context of this study on benthic-pelagic coupling, data were characterised into measures of nitrogen and carbon cycling.Nitrogen cycling parameters included the removal of bioavailable N (net denitrification rate; DNF (µmol N 2 m −2 h −1 )), the efflux of NH 4 + out of the sediment in the absence of photosynthesis (µmol N m −2 h −1 ) and the efficiency of N removal (denitrification efficiency; DE (%)) which was calculated as the percentage of DIN released as N 2 [42,43].Carbon cycling measures included C degradation rate (g C m −2 d −1 ) at the sediment surface (average of 1-3 cm depth; CD S ) calculated from the ROMA plates [37], benthic metabolism (SOC; µmol O 2 m −2 h −1 ) representing a proxy for the instantaneous rate of OM degradation [44], and a photosynthesis to respiration ratio (PR) which indicates the trophic status of the system and has been suggested to be a major control of benthic nutrient fluxes [42].PR data was taken from Mangan, et al. [41] using paired SOC and gross primary production measurements (light-dark chamber O 2 flux) within each plot, where values <1 indicate a state of net heterotrophy (more C is respired than produced) and those >1 net autotrophy (more C is produced than respired) [42].These data were previously used to highlight the differences in trophic status and thus potential C production across the 6 sites.Considering the close coupling between C supply and nutrient processing, it was necessary to use these data to understand differences in benthic nutrient fluxes.
The biochemical characteristics of DOM was assessed using three-dimensional excitation-emission matrix (3D EEM) fluorescence following Pearson, et al. [45].Four fluorescence peaks were identified in all EEMs plots (Figure S1), as reviewed by Coble [46] (see Table S1 for a description of the excitation and emission wavelengths assigned to each peak and the potential sources attributed to each).Data are presented as total humic-like and total protein-like fluorescence at sediment depths of 0-2 and 5-7 cm.Humic-like fluorescence represents a proxy for total dissolved organic carbon and degraded OM [47,48], while total protein-like fluorescence has previously been used as an indicator of water quality [49], microbial activity and has strong correlations with dissolved organic N as well as sediment O 2 consumption [50].

Statistical Analyses
While the continuation of a larger experiment by Thrush, et al. [27] allowed longer term enrichment effects to be investigated, the reduced number of sites and replication level resulted in low categorical statistical power, i.e., when separated by site and/or nutrient enrichment treatment.However, this limitation does not preclude an examination of the pooled data to test whether variability in nutrient enrichment effects was correlated to environmental variables, and whether general effects between treatments and over time could be detected.To determine if nutrient enrichment influenced sediment properties, univariate measures of macrofaunal communities and measures of N and C cycling, one-way PERMANOVAs were conducted.Nutrient enrichment (treatment) was set as a fixed factor (3 levels; control, medium and high) on data collected in the summer (n = 18).Where a winter comparison was available (n = 9), two-way repeated measures PERMANOVAs were used to investigate the effects of nutrient enrichment (fixed factor; 3 levels) and season (fixed factor; 2 levels; winter and summer), with replicate plot nested within treatment.Euclidean distance was used for all matrices.Where a significant interaction was present, the main effects were not considered and instead, post hoc tests were used to identify differences between treatments and seasons.Response parameters investigated include a multivariate analysis of sediment characteristics (organic content, mud content, median grain size), MPB biomass (chlorophyll a), phaeopigments, humiclike and protein-like DOM fluorescence at 0-2 and 5-7 cm and porewater NH 4 + and NO 3 − concentrations at 0-2 and 5-7 cm.All PERMANOVA analyses were performed using the PERMANOVA+ package in PRIMER 7.
Measures of N and C cycling were highly variable with and without the addition of nutrient enrichment.Considering the gradient of environmental variables and abundances of key bivalve species, correlations with measures of N and C cycling were explored.Pearson's correlation matrices were calculated for control, medium and high nutrient enrichment treatments separately.Variables that were depth resolved (i.e., humic-like and protein-like DOM and porewater NO 3 − and NH 4 + concentration) were averaged prior to use in the correlation analyses.
The correlation matrices revealed the relationships between N and C cycling measures and macrofaunal and environmental parameters were modified with nutrient enrichment.To explore how increasing N enrichment modified important macrofauna and environmental indices and the potential origins of these altered relationships, factor-ceiling relationships were investigated.This approach allows the assessment of any changes in response (change in the slope or shape of the ceiling) to an increasing stressor (e.g., increased nutrient enrichment) and thus implies a constraining factor where the independent variable limits the possible magnitude of the response variable [51].It also allows nonlinear relationships to be detected across all nutrient enrichment treatments.Log transformed surface (0-2 cm) porewater NH 4 + was used as a proxy for nutrient enrichment treatment, and response variables included N and C cycling and univariate macrofaunal community indices.Maximum responses at the 90th quantile were modelled as proposed by Blackburn, et al. [52].For each regression, data were divided into 4-6 bins or knots of approximately equal number and a constrained quantile curve fitted using a quadratic spline with up to two degree polynomials.The number of bins used was chosen to maximise the fit of the model, as assessed by comparing residual vs. predicted values.Data for the correlation matrices were analysed using the corrplot package and factor-ceiling relationships using the cobs package in the R software package.

Nutrient Enrichment Effects on Benthic-Pelagic Coupling
The effects of sediment nutrient enrichment on environmental properties are summarised in Table 1 (for statistical results see Table S3).Overall, nutrient enrichment and season (across a reduced number of sites; n = 3) did not significantly alter sediment properties.Univariate measures of macrofaunal communities predominately showed a negative response to nutrient enrichment unless large variability between sites (e.g., the abundance of adult A. stutchburyi ranged from 0-57 individuals per core between sites) resulted in no treatment effect, and no temporal effect was detected.Sediment nutrient enrichment significantly increased porewater N concentrations with respect to NH 4 + and NO 3 − (Tables 1 and S3).Surface (0-2 cm depth) NH 4 + was on average 40 and 800-fold higher in medium and high treatments, respectively (p < 0.001) than in control plots, which was consistent through time (p = 0.85; Table S3).Additionally, porewater NH 4 + was typically higher at depth (5-7 cm), with increases of 70 and 990-fold in medium and high treatments, respectively.Comparatively, NO 3 − concentrations were lower than NH 4 + .At the surface NO 3 − was 5 and 26 times higher in the medium and high treatments, respectively, compared to the control plots while NO 3 − at depth was 19 and 59 times higher, respectively.Inorganic nitrogen fluxes were differentially influenced by sediment nutrient enrichment (Figure 2).DNF did not differ between nutrient enrichment treatment or season (comprising a reduced number of sites (n = 3)), with rates ranging from 3.6-207 µmol N m −2 h −1 (Figure 2A, Table 2).Effluxes of NH 4 + from the sediment predictably increased with nutrient enrichment (Table 2), however, were lower during the winter compared to summer (59, 57 and 75% lower in control, medium and high treatments, respectively).Considering DNF rates remained relatively constant with enrichment and through time, reductions to the relative percentage of total DIN returned as N 2 (DE) with increasing nutrient enrichment largely reflected changes in NH 4 + efflux (Figure 2).For example, DE averaged 86, 7 and 2% in control, medium and high nutrient enrichment treatments, respectively, with DE being higher in winter compared to summer periods (Table 2).showed no treatment effects during the summer sampling (Table 2), but both were highly variable across all treatments (e.g., CDS ranged from 4-35 g C m −2 d −1 ).Parameters used to represent changes in C processing were largely unaffected by increasing nutrient enrichment (Figure 2).No difference in SOC was detected in the summer after 20 months of nutrient enrichment (Table 2).However, when comparing SOC across two seasons, a medium level of nutrient enrichment was significantly higher than both control and high treatments.Additionally, significant reductions of 62-76% were observed in winter compared to summer periods in control, medium and high treatments (Table 2).Surface carbon degradation rate (CD S ) and the trophic state of the system (PR) showed no treatment effects during the summer sampling (Table 2), but both were highly variable across all treatments (e.g., CD S ranged from 4-35 g C m −2 d −1 ).

Correlations between Biogeochemical Cycling, Environmental Variables and Macrofaunal Communities
Measures of carbon and nitrogen cycling were variable with and without the addition of nutrients, therefore correlations with environmental variables and macrofaunal communities were investigated during the summer sampling (Table 3; for full Pearson's correlation matrices see Figure S2).In control plots, DNF was most strongly correlated to C related variables (e.g., OC and H DOM |r|= 0.45-0.49),as well as porewater   The relationships and thus the contributions of particular environmental variables were altered under increasing nutrient enrichment (Table 3).In general, there was a reduction in the strength or loss of relationships between N and C cycling variables and macrofauna indices, with the only exception being an increase in strength of the correlation between DE and adult A. stutchburyi under high nutrient enrichment (r = 0.61).Instead, there was an increase in correlation strength with porewater NH 4 + concentration for the majority of N and C cycling variables (|r| = 0.45-0.69),except for CD S where no significant relationship was observed under high enrichment.
Strong interdependencies were found between N and C cycling parameters under no nutrient enrichment.In agreement with the relationships to environmental variables described above, C degradation (CD S ) was most strongly correlated to DNF and DE (r = 0.50-0.77).Nutrient regeneration, measured as NH 4 + efflux within the control treatments, was most strongly correlated to PR (r = −0.70)such that the more heterotrophic the system, the more NH 4 + was released to the water column (which in turn was strongly correlated to water column turbidity (data not shown; r = 0.64, p < 0.001)).However, no relationship between dark NH 4 + effluxes and sediment chlorophyll a was observed suggesting this was unlikely to be attributed to dark uptake by MPB.Instead, sites with lower NH 4 + effluxes were likely to have higher rates of DNF (r = −0.47).This is also revealed by the positive relationship between DE and PR, suggesting the more autotrophic the system, the higher the percentage of N was removed as N 2 (r = 0.66).
Increasing enrichment altered the relationships between N and C cycling variables.In particular, while the links between C and DNF and DE were maintained with enrichment, the relative contribution of a particular variable changed.For example, under increasing nutrient enrichment the relationship between DNF and CD s was lost but the correlation to SOC increased (r > 0.53).Additionally, the relationship between DE and CD S was weakened and then lost under high nutrient enrichment.Other changes include alterations in correlation directions, such as the reversal of relationships with NH 4 + effluxes in both medium and high treatments (Table 3).

Investigating Nonlinear Responses under Increasing N Enrichment
Considering the increasing correlation strength between N and C cycling variables to porewater NH 4 + concentration, surface (0-2 cm) porewater NH 4 + concentration was used as a measure of stress (enrichment) to examine factor ceiling relationships and thus identify any changes in response.Factor ceiling regressions between N and C cycling parameters and univariate measures of macrofaunal communities revealed a variety of different responses to enrichment, indicating specific sensitivity to increasing porewater NH 4 + concentrations (Figure 3).Response types identified included decline, increase, unimodal and skewed unimodal.NH 4 + efflux increased and DE decreased within increasing N availability in agreement to the relationships previously described (Figure 2).DNF decreased with increasing porewater NH 4 + concentration plateauing under medium and high levels of nutrient enrichment, a contrast to the comparison of DNF means (Figure 2), highlighting the loss of high variability and DNF hotspots under increasing nutrient enrichment.The skewed unimodal response of adult A. stutchburyi suggests partial tolerance to increasing NH 4 + concentrations before reaching a threshold of change resulting in linear declines.Conversely, the abundance of adult M. liliana declined with increase porewater NH 4 + such that where abundances were initially high, they reached zero at high nutrient concentrations.Factor-ceiling relationships of the trophic state (PR) of the system revealed a change in response at high levels of nutrient enrichment, where the system transitioned from being net autotrophic to increasingly more net heterotrophic, and therefore a switch from net C production to respiration.to increasing NH 4 concentrations before reaching a threshold of change resulting in linear declines.Conversely, the abundance of adult M. liliana declined with increase porewater NH 4 + such that where abundances were initially high, they reached zero at high nutrient concentrations.Factor-ceiling relationships of the trophic state (PR) of the system revealed a change in response at high levels of nutrient enrichment, where the system transitioned from being net autotrophic to increasingly more net heterotrophic, and therefore a switch from net C production to respiration.

Discussion
Overall, this study highlights how substantial alterations in the correlations of N and C cycling parameters, environmental variables and univariate macrofaunal indices under increasing nutrient enrichment can result in significant changes to the architecture of ecosystem interactions, which can significantly alter the resilience of a system to further stress [10,26,53].In addition, this study shows considerable variability in measures of N and C cycling over different temporal and spatial scales (Figure 2), which resulted in nutrient enrichment effects only being detected in a few measured end points of biogeochemical cycling (e.g., NH 4 + efflux; Table 2).However, threshold responses were observed in relation to increasing N enrichment.These included: (1) non-linear responses of macrofauna to enrichment; (2) changes in solute concentration depth profiles which may further influence macrofaunal distributions (and the associated feedbacks); and (3) photosynthesis respiration interactions.
The largest influence on biogeochemical cycling was nutrient enrichment.Enrichment of the sediment successfully increased sediment porewater concentrations of nitrogen by 40-and 800-fold in the medium and high treatments, respectively.These sediment NH 4 + concentrations are comparable to those in eutrophic estuaries globally [35], and are the same concentrations that have been used in multiple experiments investigating the role of enrichment on biogeochemical cycling, e.g., [32,38,54,55].Whilst eutrophication causes a wide variety of changes in sediment biogeochemistry including changes in C quality and quantity and O 2 concentrations, and can be dependent on differences in hydrodynamic conditions, large increases in N are commonly observed in estuaries with diffuse source inputs.For example, extensive use of urea fertilisers for vegetable cultivation around the Manukau harbour (one of our study estuaries west of Auckland city), has caused diffuse groundwater pollution, leading to nitrate-N concentrations that are ~25 mg/L [56].Many of the rivers that discharge into Manukau Harbour are predominantly groundwater-fed [57], leading to very high concentrations (~15 mg/L) of total nitrogen in-stream.
The most prominent change to nutrient processing under increasing nutrient enrichment was the increase in NH 4 + efflux.While this was reduced in the winter sampling (most likely a product of reduced activity by key bioturbators and thus the advection of solutes and microbial activity [44,58]), an overall increase in NH 4 + efflux has been linked to the beginning of a shift towards a more eutrophied state [42].This is because of potential increases in pelagic production and a reduction in N being assimilated by benthic producers [58], especially if hydrodynamic conditions encourage the retention of water column nutrients.However, DNF has the potential to mitigate or limit NH 4 + efflux owing to the removal of excess bioavailable N. DNF rates within this study were highly variable (between 4 and 208 µmol m −2 h −1 ) but comparable to those reported elsewhere in New Zealand and globally [38,42,59,60].While the factor ceiling analysis suggested a potential suppression of DNF hotspots with high porewater NH 4 + concentrations and thus the constraining of variance with enrichment (Figure 3), on average, DNF did not change with nutrient enrichment (Figure 1).Instead, higher rates of DNF within the control treatments corresponded with lower rates of NH 4 + efflux (Table 3).This is highlighted by maximum values of DNF efficiency (DE) of 220%, suggesting more N was removed as N 2 than was released back into the water column as NH 4 + (>100%).A likely consequence of these differences and thus the likelihood of NH 4 + being denitrified versus released into the water column may be associated with the trophic state of the system, such that the more heterotrophic the system the more likely NH 4 + was effluxed rather than denitrified [42].Considering nitrification (NTR) and DNF are often tightly coupled within low nutrient estuaries [61], a reduction in C production in heterotrophic systems would reduce coupled NTR-DNF through a reduction in electron donors needed for DNF, and a reduction in the supply of NH 4 + needed for NTR owing to reduced OM degradation, resulting in NH 4 + being effluxed rather than assimilated or denitrified.Further influences of nutrient stress were observed when assessing correlations between environmental variables and N and C cycling (Table 3).In particular, changes in correlation strength and direction suggests a potential switch from tightly coupled NTR-DNF to direct DNF.This would support previous evidence of a switch to direct DNF under moderate nutrient enrichment [60].Within our study, the change to direct DNF is in part a consequence of both C supply and degradation (OC, humic-like DOM and CD S ) not increasing with enrichment, suggesting DNF was likely to be C limited.Additional indicators of this potential change include increases in the correlation strength of DNF to porewater nutrient concentration, and reductions in strength to chlorophyll a and macrofaunal abundances with increasing enrichment.The facilitation of direct DNF would therefore manifest as a consequence of the increase in nitrate concentration reducing the competition between primary producers and nitrifying and denitrifying bacteria [62].
Alterations in the correlation strength and direction of N and C cycling, environmental variables and macrofauna indices were driven by threshold responses in key relationships as N stress increased (Figure 3).For example, the unimodal response of a dominant macrofauna species, A. stutchburyi, to increasing N enrichment demonstrates a change in response and reductions in total abundance under a medium level of nutrient enrichment, highlighting sensitivity to increasing N enrichment.The sensitivity of M. liliana was significantly greater than A. stutchburyi, with linear declines observed with increasing enrichment.These relationships to nutrient enrichment are typical of those reported elsewhere [23,24,63].While initial increases in adult A. stutchburyi may facilitate coupled NTR-DNF through increased bioturbation activity [64], eventually, reductions in the abundance of both key bioturbating species will reduce the spatial and temporal heterogeneity of sedimentary redox zones and thus suppress NTR under increased sediment anoxia [64][65][66][67].
The reduction in the vertical movement of solutes within the sediment with increasing nutrient enrichment was evidenced by comparing the difference between surface and deep porewater DOM under different nutrient enrichment treatments (Figure S1).Humic-like DOM at the sediment surface was 6% lower than at depth in both control and medium nutrient enrichment treatments but 33% lower in response to high enrichment.The larger difference at high N enrichment indicates a more stratified sediment layer and thus lower mixing and bioturbation activity.Furthermore, the higher intensity of humic-like DOM at depth has been suggested to represent the accumulation of refractory, low molecular weight DOM during sediment OM diagenesis in addition to increased sediment anoxia [47].These increases sediment anoxia may have further contributed to the reductions in macrofauna abundances and thus bioturbation activity in the high nutrient enrichment treatments.
A unimodal relationship was additionally observed for the trophic status of the system (as indicated using PR), where a critical transition in response occurred between medium and high nutrient enrichment (Figure 3).The modification of PR is an important indicator of benthic nutrient fluxes in shallow coastal ecosystems because increasing net heterotrophy can lead to an increase in the release of inorganic nutrients through respiration [42].Furthermore, the preceding study by Thrush, et al. [27] which evaluated changes in ecosystem interaction networks associated with nutrient processing for the initial 9 month enrichment period, reported that high water column turbidity (and thus heterotrophy) resulted in a simplified ecosystem interaction network that included the decoupling of sediment nitrogen processing and thus a reduced capacity to process excess nutrients.
The transition to increasing net heterotrophy has important implications for overall ecosystem functioning owing to the potential intensification of pelagic production and the start of a cascade towards a more eutrophic state.These non-linear relationships are particularly important considering a small increase in cumulative stress can result in an abrupt ecological shift if thresholds are crossed [68].
It is important to recognise that the patterns presented within this study are correlative and there are other factors which can influence sediment biogeochemical cycling.For example, concomitant increases in sediment phosphorus (P), which were not stimulated in this enrichment study, may have altered nutrient dynamics through an increased stimulation of primary producers [69] which in turn influences sediment oxygen concentrations and the recycling and regeneration of N within the sediment.While P within this study was close to detection limits, it is likely that an increase in the release of P from the sediment may have been observed under combined N and P enrichment, especially under increasing sediment anoxia [70].Other factors influencing biogeochemical cycling include changes in quality and quantity of any accompanying C inputs [71] in addition to the stimulation of hypoxia which often occurs simultaneously with coastal eutrophication [22].Nonetheless, the importance of understanding the response of multiple estuaries to increasing N loading helps to further our understanding of the changes to structure and functioning of coastal communities to point source N pollution.

Conclusions
Both C and N cycling within oligotrophic estuarine sediments are tightly coupled and embedded within a complex network of ecosystem interactions, involving direct and indirect feedbacks.However, anthropogenic pressure on coastal ecosystems is intensifying which is impeding the capacity to moderate intensifying pollutant loads.By coupling measurements of N and C cycling, this study highlights the limited capacity for nitrogen removal in response to large increases in N availability and the loss of DNF hotspots under increasing enrichment, as evidenced through the reductions in DE and a reduction of DNF with increasing porewater NH 4 + concentration.Increases in porewater N concentrations along with subsequent changes in macrofaunal communities resulted in the potential uncoupling of NTR-DNF to direct DNF.Macrofaunal communities and the degree of net autotrophy/heterotrophy (PR) were shown to have important influences on N processing.Factor-ceiling relationships revealed, that in general, maximal abundances of macrofauna corresponded to lower N enrichment, highlighting the vulnerability of estuaries to increasing nutrient loads.The appearance of these transitions can signify changes to the interactions of intrinsic dynamics and drivers which can fundamentally alter biogeochemical cycling within soft sediments and increase the likelihood of abrupt non-linear shifts in ecosystem functioning [27,72].Consequently, current limit setting management approaches to increasing N enrichment may not be suitable in preventing a shift towards a more eutrophied state as dynamic changes dominate ecosystem interactions.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jmse10121955/s1.Figure S1: 3D excitation-emission matrices (EEM) of porewater DOM at 0-2 cm and 5-7 cm sediment depth as a function of nutrient enrichment treatment.Table S1: 3D EEMs fluorescence peak descriptions.Table S2: Sediment properties and univariate indices of macrofaunal communities as a function of sediment nutrient enrichment and season.Table S3: Results of one-way PERMANOVAs to determine treatment effects on sediment and macrofauna properties, DOM fluorescence and porewater nutrient concentrations.

Figure 1 .
Figure 1.(A,B) Locations of the 6 sites, within 4 estuaries situated within the North Island of New Zealand.(C) Schematic of the plot layout within each site and (D) dates of sampling and nitrogen additions carried out after all sampling was completed.

Figure 1 .
Figure 1.(A,B) Locations of the 6 sites, within 4 estuaries situated within the North Island of New Zealand.(C) Schematic of the plot layout within each site and (D) dates of sampling and nitrogen additions carried out after all sampling was completed.

Figure 2 .
Figure 2. Inorganic nitrogen (A,C,E) and carbon related (B,D,F) fluxes and indices during summer (yellow, n = 18) and winter (blue, n = 9) as a function of nutrient enrichment treatment.Boxes represent the 25th and 75th percentiles, and whiskers are the 5th and 95th percentiles.A solid line within each box denotes the median.Note the log10 scale of the y-axis in panel C. DNF = net denitrification rate; NH4 + = NH4 + efflux; DE = denitrification efficiency; SOC = sediment oxygen consumption; CDS = carbon degradation rate at the sediment surface; PR = photosynthesis respiration ratio.

Figure 2 .
Figure 2. Inorganic nitrogen (A,C,E) and carbon related (B,D,F) fluxes and indices during summer (yellow, n = 18) and winter (blue, n = 9) as a function of nutrient enrichment treatment.Boxes represent the 25th and 75th percentiles, and whiskers are the 5th and 95th percentiles.A solid line within each box denotes the median.Note the log 10 scale of the y-axis in panel C. DNF = net denitrification rate; NH 4 + = NH 4 + efflux; DE = denitrification efficiency; SOC = sediment oxygen consumption; CD S = carbon degradation rate at the sediment surface; PR = photosynthesis respiration ratio.

Figure 3 .
Figure 3. Bivariate scatter plots of nitrogen (A,D,G) and carbon (B,E,H) cycling parameters and univariate macrofaunal indices (C,F,I) showing factor-ceiling relationships (black line) at the 90th percentile.Where the line is missing, the relationship did not change with increasing enrichment.Blue, grey and orange dots are given as a reference to control, medium and high nutrient enrichment treatments, respectively.Data included are from summer only (n = 54).

Figure 3 .
Figure 3. Bivariate scatter plots of nitrogen (A,D,G) and carbon (B,E,H) cycling parameters and univariate macrofaunal indices (C,F,I) showing factor-ceiling relationships (black line) at the 90th percentile.Where the line is missing, the relationship did not change with increasing enrichment.Blue, grey and orange dots are given as a reference to control, medium and high nutrient enrichment treatments, respectively.Data included are from summer only (n = 54).
Figure S2: Pearson's correlation coefficients between summer N and C cycling parameters, environmental variables and univariate indices of macrofaunal communities.

Table 1 .
Summary of sediment properties and univariate indices of macrofaunal communities as a function of sediment nutrient enrichment treatment.Data presented are the median (the range is in brackets) across winter and summer (n = 27), except for macrofaunal abundance, richness and diversity where only summer data were available (n = 18).

Table 2 .
Results of one-way PERMANOVAs comparing summer univariate measures of nitrogen and carbon processing as a function of nutrient enrichment treatment (n = 18).Repeated measures two-way PERMANOVAs are used to compare summer and winter values at three sites (n = 9).Significant effects (p < 0.05) are given in bold and post hoc pairwise tests are shown for significant interactions.Nutrient enrichment treatments have been abbreviated to C = control, M = medium, H = high, and winter and summer to W and S, respectively.

Table 3 .
Pearson correlation coefficients between (A) measures of N and C cycling and environmental variables and (B) N and C cycling during the summer as a function of nutrient enrichment treatment (n = 18).Data shown are significant relationships (p < 0.05), see Figure S2 for all relationships.* p < 0.05; ** p < 0.01; *** p < 0.001.