Water Column Turbidity Not Sediment Nutrient Enrichment Moderates Microphytobenthic Primary Production

Soft sediment intertidal habitats are under intense anthropogenic pressure resulting from increased land derived sediment and nutrient delivery. Long term, this can cause high water column turbidity and nutrient enrichment of sediment porewaters, which has cascading effects on coastal ecosystem functionality. However, how these stressors may interact and influence benthic productivity over alternating periods of submergence and emergence is largely unknown. This study investigates the effects of sediment nutrient enrichment (at three levels for 20 months) on benthic primary production at six sites in four New Zealand estuaries that spanned a gradient in water column turbidity. While nutrient enrichment had no detectable effect on microphytobenthic primary production, water column turbidity had a significant influence, explaining up to 40% of variability during tidal submergence, followed by temperature and sediment characteristics. In addition, negative net primary production (NPP) estimates and therefore net heterotrophy for the most turbid estuaries during tidal submergence resulted in an increased reliance on production during emerged periods, where NPP was positive across all sites. This study highlights the prominent role of water column turbidity over nutrient enrichment in moderating microphytobenthic production, and the increasing importance of emerged periods to maintain the health and functioning of coastal habitats.


Introduction
Soft sediment intertidal habitats are under intense anthropogenic pressure which is diminishing ecological functioning and the delivery of ecosystem services upon which humanity relies. Globally, the most pervasive pressures include significant increases in nutrient and sediment delivery [1] predominately through discharges of surface run-off and groundwater inputs [2,3]. These act to modify coastal nitrogen cycling as well as increase light attenuation to the benthos [4,5]. Ultimately, this can diminish benthic primary productivity which cascades to alter global carbon budgets, food web dynamics and water quality parameters [6][7][8]. Despite its importance, there is currently a limited understanding of how these pervasive pressures may interact and influence benthic productivity on intertidal flats over alternating periods of high and low tide within ecologically relevant contexts. Here for the first time, we couple in situ submerged and emerged primary production estimates across two seasons, in response to elevated sediment nutrient enrichment, over a natural gradient of water column turbidity.
Microphytobenthos (MPB) are often the dominant primary producer in shallow temperate ecosystems, accounting for up to 50% of total estuarine autochthonous primary production and up to 80% of total benthic carbon fixation [9]. While MPB productivity can supply high-quality labile carbon, and thus underpin coastal food webs [6], MPB also modify sediment stability and, through oxygenation of the sediment and nutrient uptake, alter biogeochemical pathways including nutrient recycling [8,10,11]. MPB are therefore fundamental constituents of intertidal habitats, with reductions in productivity likely to have cascading implications for entire coastal ecosystems.
Rates of MPB production are moderated by a number of factors such as hydrodynamic conditions, carbon supply, the availability of nutrients and light climate [12,13]. Of particular importance during submerged tidal periods is light availability [14]. Increases in water column light attenuation are predominately a product of suspended sediment concentrations [15] but can additionally be influenced by phytoplankton blooms in response to nutrient enrichment [16]. High water column turbidity restricts the amount of light reaching the seafloor and thus can directly diminish MPB production. While MPB have been shown to adapt to light conditions as low as 2.8 µmol m −1 s −1 [17], maximal rates of productivity are often proportional to light saturation values [14] and therefore productivity and the ecosystem services provided by MPB are often reduced with decreasing light availability. For example, reduced productivity (when light attenuation is high) causes feedbacks within the system via a reduction in labile organic carbon available to consumers, increases in the efflux of dissolved inorganic nitrogen from the sediment [18], reductions in macrofaunal diversity [19] and alterations to sediment stability [11].
Within subtidal environments where turbidity can remain consistently high, a tipping point can occur where the system moves from benthic to pelagic dominated primary production [20]. However, for intertidal environments, water column turbidity is a temporally displaced stressor as the tide recedes and exposes the sediment surface to high light availability. Therefore, exposure associated with low tide periods can provide some resilience to reduced MPB production during submerged periods. This has been demonstrated in several estuaries throughout the world, where benthic primary productivity is often reported to be limited to low tide periods [21][22][23][24]. The dependence on this exposure period and the relative contribution of low tide MPB production may therefore be closely coupled to site turbidity [25]. Despite this crucial link, there are few studies where direct comparisons of emerged and submerged primary production exist [23,26,27], especially through time and along a transitional gradient of water column turbidity. These knowledge gaps are limiting our ability to understand the effects of increased anthropogenic pressure.
Increased water column turbidity often occurs concomitantly with increased nutrient inputs. While it is well known that increased nutrient availability can lead to phytoplankton blooms, increases in filamentous benthic macroalgae and ultimately coastal eutrophication [28,29], the response of MPB production is less understood. This may be partially attributed to the majority of research focusing on the individual influences of elevated nutrients during one phase of a tidal cycle. However, nutrient responses by MPB are likely to be tightly coupled to light availability [30] and differ depending on tidal state. For example, an increase in nutrient availability has been shown to stimulate MPB productivity in a nutrient-limited system with high light availability during tidal submergence [31]. However, in a highly turbid but nutrient-rich estuary, light limitation was sufficient to prohibit MPB production and thus nitrogen uptake [32]. Conversely, during tidal emergence in a high nutrient system, MPB were not observed to be nutrient-limited [33,34], while in a low-nutrient system, the response of MPB is still largely unknown. Considering the multifaceted response of MPB within intertidal habitats to nutrient enrichment, it is vital for investigations to include a comparison of responses across tidal cycles, in combination with other pressures (such as water column turbidity) and which encompass timeframes that allow the incorporation of longer-term responses (e.g., microbial and macrofaunal community responses).
The objective of this study was to investigate the interactive effects of sediment nutrient enrichment and water column turbidity on both submerged and emerged MPB primary production. This was carried out by experimentally enriching the sediment at three levels for up to 20 months along a gradient of water column turbidity. Within the main objective, we postulate two questions; (1) how does water column turbidity influence submerged and emerged MPB primary production, and is this modified by sediment nutrient enrichment? In addition, (2) are these responses temporally variable? This research aims to fill a critical gap in our understanding of how two pervasive stressors may interact within natural environments, with a focus on MPB production, an often forgotten but integral component of coastal ecosystems [35].

Study Sites
This study was carried out at 6 sites within 4 estuaries of the North Island of New Zealand (Figure 1; for site-specific GPS coordinates see Table A1). All estuaries were shallow, barrier enclosed coastal systems with extensive intertidal area and had semi-diurnal tides. Sites were chosen based on a perceived gradient of water column turbidity while having a similar latitude to normalise for daylength and temperature. In addition, sites were located within the mid-intertidal region and inhabited by the functionally important bivalve species Austrovenus stutchburyi and Macomona liliana, which have been shown to significantly influence benthic-pelagic coupling [36].  Table A1.

Experimental Design
The experimental period spanned a total of 20 months with sampling in November 2017 (T1), June 2018 (T2; 3 sites only-TAU-T, TAU-O, RAG) and November 2018 (T3), corresponding to 8, 15 and 20 months of sediment nutrient enrichment (Table 1). Two November sampling periods were chosen to correspond with warmer late spring temperatures and higher light availability, and allowed the incorporation of potentially longer-term responses to nutrient enrichment (e.g., alterations of macrofaunal communities). An additional winter sampling (June 2018) at three sites was conducted to examine the potential influence of reduced temperatures and light conditions on primary production measures. At each site, nine 9 m 2 experimental plots were divided into three nutrient enrichment groups, either control 0 g N m −2 , medium 150 g N m −2 or high nutrient treatment 600 g N m −2 ( Figure 1F). To achieve enrichment, a slow-release urea fertiliser (Nutricote 40:0:0 N:P:K) was added at two known quantities to a series of 180 evenly spaced holes at a depth of 15 cm by extracting a sediment plug (3 cm diameter) with a handheld corer, adding a known volume of fertiliser and then replacing the sediment plug. Nitrogen-only fertiliser was used since New Zealand estuaries are typically N limited and because the urea in Nutricote quickly hydrolyses to ammonium (NH 4 + ), a form of nitrogen associated with eutrophication and the remineralisation of organic matter [37]. This method results in an even elevation of porewater NH 4 + throughout the experimental plots [38]. To maintain enrichment throughout the experimental period, fertiliser was first applied in March 2017 and reapplied following sampling in November 2017 and June 2018 (Table 1). Sampling, therefore, took place between 5 and 8 months after each application of fertiliser.

Primary Production Measurements
Primary production was measured during submergence at all sites in November 2017 and 2018, and at TAU-O, TAU-T and RAG in June 2018. On each sampling date, an Odyssey PAR (photosynthetically active radiation) logger (Dataflow Systems, Melbourne, FL, USA) was deployed in the middle of each site and on the shore to capture light availability at the sediment and water surface, respectively. Two benthic flux chambers were deployed in each plot for approximately 4 h during a midday high tide. Each chamber base (50 cm × 50 cm × 15 cm) was inserted 5 cm into the sediment at low tide and equipped with a light and temperature logger (HOBO Pendant ® ). On the incoming tide, chambers were sealed with a Perspex dome lid and any air bubbles removed, encapsulating 41 L of seawater over the sediment surface. An opaque shade cloth covered one of the two chamber lids per plot to block out all sunlight, while the other lid was left uncovered to receive ambient sunlight. Immediately after the chambers were sealed and at the end of the incubation, one 60 mL seawater sample was collected from each chamber to provide initial and final dissolved oxygen concentrations. Seawater dissolved oxygen concentration was measured on-site immediately after sampling using a handheld optical probe (PreSens FIBOX 3 LCD trace v7). In June and November 2018, this was supplemented with a dissolved oxygen logger (PME miniDOT) deployed within each chamber recording at a 1-min sampling frequency.
Emerged primary production was assessed within 10 d of the submerged measurements in June 2018 at TAU-O, TAU-T and RAG, and at all sites in November 2018 over a midday low tide. For the full emerged period, an Odyssey PAR light logger was deployed in the middle of the site to capture ambient light availability. Once the site was drained of water, one chamber base in each plot was inserted 5 cm into the sediment. During each measurement, the chamber contained an Odyssey PAR light logger, a pressure vent, a thermocouple measuring chamber air temperature and a battery-powered fan. The fan maintained airflow and gentle mixing to ensure no dead spaces were created inside the chamber which could alter the diffusion of gases from the sediment [39]. A Perspex lid was then fitted over the chamber and connected to an infrared gas analyser (LI-COR 8100A Automated Soil CO 2 Flux System) where air was continuously circulated between the analyser and the chamber. Measurements of CO 2 concentration (ppm) were recorded at a frequency of 1 Hz for a total incubation period of 5 min. This method has been shown to provide a reliable quantification of CO 2 exchange in intertidal habitats, with a 5 min period resulting in a stable diffusion of CO 2 without increasing the humidity within the chamber [25,40]. Light incubations were conducted before an opaque shade cloth was placed over the chamber base and allowed to acclimate for 20 min prior to an additional CO 2 incubation performed in dark conditions.

Site Characteristics
Water column turbidity was assessed by measuring light availability reaching the sediment surface every 10 min over a 9-month period (between March and November 2017; see Mangan, et al. [14]), with PAR Site an average of the maximum daily submerged PAR values during this period. To characterise sediment properties, on each submerged sample date, five pooled surface sediment cores (2.6 cm dia, 2 cm deep) were taken from each experimental plot and frozen at −20 • C [41] until analysis for sediment chlorophyll a, phaeopigments, grain size and organic matter content. Two separate replicates of 4 additional sediment core samples at depths of 0-2 cm and 5-7 cm were taken within each plot and kept on ice for analysis of sediment porosity and porewater ammonium concentrations. To assess any potential desiccation during emerged primary production measurements, four pooled 0-2 cm and 5-7 cm sediment cores (2.6 cm dia) were also taken and kept on ice until analysis of sediment porosity. Two macrofauna cores (13 cm dia, 15 cm depth) were taken from each plot in November 2017 and 2018 and sieved on a 500 µm mesh before being preserved in 70% Isopropyl alcohol.

Laboratory Analysis
Sediment samples were thawed and homogenised before analysis. To determine concentrations of chlorophyll a and phaeopigments, a sub-sample of sediment was freeze-dried, extracted with 90% buffered acetone and measured before and after acidification using a fluorometer (Turner Designs 10-AU). Grain size was measured using laser diffraction (Malvern Mastersizer-3000) after the digestion of organic matter with 10% hydrogen peroxide. Organic content was determined by weight loss on ignition following 3 d at 60 • C to ensure a constant weight and then combustion at 550 • C for 4 h. Sediment samples collected for porewater ammonium determination underwent water extraction within 24 h of collection. Four ml of de-ionised water was added to each sample before being vortexed and centrifuged. The supernatant was then filtered through a 0.45 µm Whatman GF/C glass fibre filter and stored at −20 • C. Analysis for porewater ammonium concentration was conducted on a Lachat QuickChem 8000 Series FIA+ (Zellweger Analytics Inc. Milwaukee, Wisconsin 53218, USA) using standard operating procedures for flow injection analysis. Sediment porosity samples (both for porewater calculations and those taken during each emerged incubation) were analysed by calculating the difference in wet and dry weight after 7 days at 60 • C (or until constant weight). Macrofauna samples (6 per treatment in T1 was reduced to 3 per treatment in T3 as no statistical difference in univariate measures were detected) were stained with Rose Bengal and fauna separated from any sediment and shell hash before being identified to the lowest possible taxonomic level (usually species) and counted.

Data Analysis
Submerged dissolved O 2 and emerged CO 2 fluxes were estimated from the change in concentration during the incubation period and corrected for the chamber area and volume. Submerged O 2 fluxes were not corrected for water column processes owing to the small contribution relative to benthic production/respiration (<5% of the measured chamber flux, data not presented). A respiratory quotient of 1 was used to convert CO 2 consumption and production measurements to O 2 production and consumption, respectively. This value is likely to be realistic for well-oxygenated sediments with low mud and organic content [42], such as those within this study. Light chamber fluxes were equivalent to net primary production (NPP; µmol O 2 m −2 h −1 ) and dark chamber fluxes to sediment oxygen consumption (SOC). The sum of these two fluxes from paired chambers within a plot provided an estimate of gross primary production (GPP). A productivity/respiration (p/r) ratio was calculated for each chamber, as defined by Eyre and Ferguson [43]: (hourly productivity (GPP) × daylight period)/(hourly respiration × 24 h), to estimate if the sediments were net autotrophic (p/r > 1; more carbon is produced than respired), or net heterotrophic (p/r < 1; more carbon is respired than produced).
Odyssey PAR logger data were converted from count measurements to PAR (µmol m −2 s −1 ) by a calibrated regression with a LI-COR PAR sensor prior to deployment. Three values of PAR are used throughout this study: site PAR (PAR Site ), incident PAR (PAR i ) and sediment surface PAR (PAR Sed ). PAR Site is used as a proxy for site turbidity as described above, where a higher PAR Site value indicates lower site turbidity. The other two PAR measurements correspond to those taken during sampling and are used for normalising rates of primary production. PAR i describes light received at the sediment surface during emerged incubations and at the water surface during submerged incubations. Whereas PAR Sed refers to light received at the sediment surface during the submerged incubations and therefore accounts for water column attenuation. GPP was first normalised by PAR i (GPP i ) to account for variability in light intensities between sampling dates while retaining attenuation effects of the water column (e.g., from turbidity). GPP was additionally normalised by PAR Sed and photosynthetic biomass (sediment chl a i.e., productivity per unit of photosynthesising biomass [18]) to determine photosynthetic efficiency (GPP Sed + chla ).
To determine if nutrient enrichment influenced primary production estimates and if this differed between sites and sampling events, a preliminary analysis using repeated measures PERMANOVA was conducted. For submerged fluxes, nutrient enrichment (3 levels), site (6 levels) and time (2 levels; T1 and T3) were set as fixed factors and replicate plot nested within treatment. T2 was omitted because measurements were only made at 3 of the 6 sites. For emerged fluxes, site was reduced to 3 levels (RAG, TAU-O, TAU-T) and time included T2 and T3 ( Figure A1, Table A5). These analyses revealed high spatial and temporal variability. This was expected due to the extent of the experiment (6 different sites within 4 estuaries, over 20 months) which spanned a range of environmental conditions expected to influence primary production. To differentiate between this natural heterogeneity (in both time and space) and reveal any potential treatment effects, submerged and emerged primary production estimates (NPP and GPP) were normalised by an average of the control plot values. This resulted in treatment effect size being relative to the site-specific background level and sampling date. Control normalised (CN) submerged and emerged primary production estimates were tested for differences from control values (i.e., NPP CN , GPP CN 1) using one-sample t-tests. Differences between nutrient enrichment treatments (medium vs. high) were tested using two-sample t-tests and between sampling periods (for submerged data only) (T1 vs. T3) using paired two-sample t-tests. Statistical analyses for the control normalised data were performed using the R stats package (v3.6.2) in R Studio.
No significant treatment effects on primary production were detected in the PERMANOVA ( Figure A1, Table A5) or control normalised analyses (see Results) and therefore raw (i.e., not control normalised) data were pooled for subsequent analysis. To investigate if site turbidity modified the relative importance of submerged and emerged primary production and to see if this was consistent through time, two repeated measures PERMANOVA's were conducted. The first using T3 (November 2018) data with site and tidal state as fixed factors (6 and 2 levels, respectively) and replicate plot nested within site. The second considered a temporal element by including site (fixed factor, 3 levels: TAU-O, TAU-T and RAG), time (fixed factor, 2 levels: T2 and T3) and tidal state (fixed factor, 2 levels: submerged and emerged) with replicate plots nested within the site. All similarity matrices were based on Euclidian distances, and primary production measures tested were NPP, GPP, GPP i , GPP Sed + chla . Main effects were not considered if interaction terms were significant, and instead, posthoc tests identified differences between sites, tidal state and sampling date.
Distance-based Linear Models (DistLM) were then used to identify if any environmental (sediment characteristics, light etc.) or univariate macrofaunal variables were significant drivers of primary production measures across all nutrient enrichment treatments, sites and sampling dates. First, predictor variables were normalised to allow for comparisons despite differing units and scales. Significant individual predictors (marginal tests) could then be identified and the best combination of predictor variables (stepwise procedure) of NPP and GPP during submerged and emerged tidal periods calculated. The corrected Akaike information criterion (AICc) was used for all models as AICc is suggested to be most appropriate when there are a large number of predictor variables relative to the sample size (sample number/n. explanatory variable < 40) [44]. Where co-linearity occurred among predictor variables (r > 0.7), the variable explaining the least amount of variability was omitted from the full model. For June 2018 (T2), macrofauna data between the two November samples were averaged so this sample date could be included. Variance partitioning was used to determine how much of the observed variability was attributed to significant environment variables (organic content, mud content, median grain size, sediment chl a content, phaeophytin content, PAR i , PAR sed (submerged data only), PAR site and temperature (sediment surface)), community variables (total abundance, taxa richness, A. stutchburyi abundance and M. liliana abundance) and porewater variables (NH 4 + concentration at 0-2 cm and 5-7 cm sediment depth and porosity at 0-2 cm (emerged data only)). The two bivalve species were included based on their known influence on benthic processes within these systems [36]. All PERMANOVA and DistLM analysis were performed using PRIMER 7 with the PERMANOVA+ package.

Site Characteristics
During the deployment of the PAR light loggers, no fouling or sediment deposition was observed and therefore the difference in measured light availability between sites during submergence confirms a gradient in water column turbidity was captured across the six sites (PAR Site ; Figure 2A). Light availability was lowest at MNK-L and MNK-R with a median maximum daily PAR of 218 and 290 µmol m −2 s −1 , respectively. In addition, the interquartile range (IQR) was relatively low for both sites (297 and 365 µmol m −2 s −1 , respectively), which is suggestive of a consistently low light climate and thus high turbidity. RAG while having a similar IQR to both the MNK sites (283 µmol m −2 s −1 ), had a higher median maximum daily PAR of 500 µmol m −2 s −1 . In contrast, TAU-O, TAU-T and WTA had a considerably higher median maximum daily PAR (784-1140 m −2 s −1 ), suggestive of overall less turbid water, and a high IQR (512-762 m −2 s −1 ), indicative of intermittent increases in turbidity (rather than persistently high). Sediment surface PAR (PAR Sed ) measurements during submerged sampling periods suggest this gradient of water column turbidity was captured during each sampling event ( Figure 2B-D). Where high incident light (PAR i ) occurred, greatest differences between PAR Sed were observed at MNK-L, MNK-R and RAG. However, these differences are less pronounced with lower PAR i values (e.g., during T2 sampling in winter). No consistent relationship was found between water column turbidity and sediment characteristics ( Table 2). On average, mud content was relatively low (<4%) at MNK-R, TAU-T and WTA on all sampling dates. The highest mud content was recorded at RAG (maximum of 28%) with both MNK-L and TAU-O predominately >10%. Organic content was also relatively low at all sites (<5%) while median grain size (MGS) varied between sites (an average of 123 µm in RAG to 248 µm in WTA). Microalgal biomass (sediment chl a content) was highly variable with no clear relationship between sites or sampling dates. Univariate measures of macrofaunal community differed between sites ( Table 2; see Tables A3 and A4 for community composition data). Most notable was the high abundance of A. stutchburyi (>50 core −1 ) at MNK-L, RAG and WTA which contributed to higher total macrofaunal abundance.   OC = organic content, Mud = mud content, MGS = median grain size, Chl a = chlorophyll a, Phaeo = phaeopigment, DW = dry weight; N = total abundance, S = taxa richness.

Treatment Effects on Sediment Properties
Sediment nutrient enrichment successfully increased porewater ammonium concentration at all sites throughout the experimental period (Table A2). On average, porewater ammonium (NH 4 + ) increased 40-fold in medium treatments and 800-fold in high treatments relative to control concentrations at a depth of 0-2 cm. Porewater NH 4 + concentration was typically higher at a depth of 5-7 cm relative to 0-2 cm, with 70 and 990-fold average increases in the medium and high treatments, respectively (relative to control concentrations). Successful nutrient enrichment consequently altered sediment porewater N:P ratios ( Figure A2). Within control plots, N-limitation was predicted to occur at all sites through time (N:P < 16), while a switch to P-limitation likely occurred in all enriched plots through time (assuming a switch from N-to P-limitation is likely to occur at an N:P ratio between 10 and 16 [45,46]). Overall, sediment enrichment did not translate into substantial changes in sediment properties (grain size, mud content and organic content) (Table A2). A medium level of nutrient enrichment typically did not alter macrofaunal communities (Tables A2-A4). However, negative responses were observed within high treatment plots, such as reduced abundance or loss of A. stutchburyi and M. liliana (Table A2).

Treatment Effects on MPB Biomass and Primary Production
MPB biomass varied with site, sample date and sediment enrichment level, but not in a consistent manner (e.g., high treatments having greater biomass than other treatments, or consistent relationships within a site through time) (Table A2). Any deviations in MPB biomass from control plots primarily occurred at sites with the lowest water column turbidity (TAU-O, TAU-T and WTA) while no difference in chl a was observed at the most turbid sites (MNK-L and MNK-R during T1 and RAG throughout the experimental period) (Table A2).
Preliminary analysis of nutrient enrichment effects on primary production measures revealed high variability between sites and sample dates ( Figure A1, Table A5), and so site-specific data were control normalised to see whether this would reduce variability and reveal effects of nutrient enrichment on primary production (Figure 3). Nutrient enrichment did not have a significant effect on submerged NPP CN or GPP CN after 8 (T1) and 20 (T3) months of enrichment regardless of treatment level (medium or high) ( Table 3). While it was not possible to distinguish long-term treatment effects (between T1 and T3) for emerged production measures, nutrient enrichment effects on NPP CN and GPP CN were not observed after 20 months of enrichment (T3). Therefore, nutrient enrichment treatment data were pooled for all subsequent analysis.

Submerged vs. Emerged Primary Production and Its Predictors
To investigate the relative importance of submerged and emerged primary production and if this differed over a gradient of water column turbidity, primary production measures as a function of site and tidal period during T3 were compared (Figure 4). At the three sites with the highest water column turbidity (MNK-L, MNK-R and RAG), submerged NPP was significantly lower than emerged NPP and negative at the two most turbid sites (MNK; Table 4). In contrast, at the three sites receiving the greatest light availability, submerged NPP was equal to or greater than emerged NPP. The relationship between GPP and water column turbidity was weaker, with emerged GPP at the three most turbid sites being only slightly greater than or equal to emerged GPP from the least turbid sites ( Figure 4B). After standardising GPP by incident PAR, submerged GPP i did not strongly relate to water column turbidity. This may be partially attributed to the low incident light recorded at RAG ( Figure 2D). Emerged GPP i was highly variable among sites and did not significantly exceed submerged GPP i (Table 4). Emerged GPP i was considerably higher at TAU-O when compared to the other sites, however, this is likely to be an artefact of both low light availability during the emerged sampling, and high incident light during the submerged sampling event ( Figure 2D). Photosynthetic efficiency (GPP Sed + chla : correcting GPP by PAR Sed and photosynthesising biomass (chl a)) appeared to be independent of water column turbidity during submerged tidal periods, with little difference observed among sites ( Figure 4D). The only exception occurred at MNK-R where photosynthetic efficiency was notably higher. Submerged GPP Sed + chla was greater than emerged GPP Sed + chla at all but one site (TAU-O) where photosynthetic efficiency was equal between tidal periods. These relationships appear to be relatively robust through time at the three sites where a June comparison is available ( Figure 5). Submerged NPP was lowest at the most turbid site (RAG) during T2 and T3, and significantly lower than emerged NPP compared to sites with higher light availability (TAU-O and TAU-T). Emerged NPP at RAG was significantly higher than both TAU sites during both sampling dates (Table 5). All GPP estimates showed no clear relationship with time, however, RAG was observed to have the lowest submerged GPP and highest emerged GPP across both sampling dates. Once corrected for incident light, submerged GPP i was three times higher at TAU-O, compared to all other sites and across the two sampling dates, while all other submerged GPP i fluxes were relatively consistent. Photosynthetic efficiency remained significantly higher during submergence at all sites and sampling dates except for RAG during T2 and TAU-O during T3 where similar estimates between tidal periods were observed.    The influence of water column turbidity on overall productivity at each site through time is highlighted in Figure 6. Net autotrophy (more carbon produced than respired) dominated all sites during tidal emergence. However, during submerged tidal periods, sites with high water column turbidity consistently had a p/r ratio <1 suggesting net heterotrophy (more carbon respired than produced). In contrast, sites with low water column turbidity predominately remained net autotrophic during submergence at all three sampling dates.
The variability in primary production among sites and tidal stages was partitioned using DistLM analysis ( Table 6). The analysis revealed environmental variables contributed 10-25% of overall variability for submerged and emerged primary production estimates. When considering the shared effects of community variables for submerged estimates, and porewater variables for emerged estimates, the total amount of explained variability increased to 19-39%. PAR Site was significantly correlated with submerged primary production estimates, explaining 8-40% of the total variance when considered individually, supporting earlier evidence of the influence of water column turbidity on benthic primary production. Temperature explained the second largest proportion of variability and was included in the full model during both submerged GPP and emerged NPP estimates (11% and 4%, respectively). Other significant individual predictors for both submerged and emerged estimates included sediment properties (either mud content, organic content or grain size), phaeopigments and porewater NH 4 + concentration while the addition of total abundance of A. stutchburyi and porosity were significant for emerged estimates. Porewater NH 4 + concentration explained a maximum of 4% variability, further highlighting the marginal influence of sediment enrichment on primary production estimates. Figure 6. Photosynthesis/respiration ratio (p/r) during submerged (blue) and emerged (red) tidal periods with data pooled across all sampling dates and treatments. Boxes are the 25th and 75th percentiles, and whiskers the 5th and 95th percentiles. The dotted line is provided as a reference to p/r = 1.

Discussion
Comparisons of MPB productivity on intertidal flats during periods of tidal inundation and exposure are rare, despite the potential of several anthropogenic stressors to restrict productivity to emerged periods. This study considered both submerged and emerged productivity and integrated these estimates over a 20-month period within six different estuaries, to investigate the combined effects of two pervasive coastal pressures-increased nutrient availability and water column turbidity. To achieve this, MPB productivity was measured at both high and low tide (i.e., photosynthetic O 2 production and CO 2 uptake on submerged and exposed tidal flats, respectively) and analysed according to nutrient enrichment treatment across a natural turbidity gradient (Figure 2A). Experimental nutrient enrichment resulted in sediment porewater ammonium concentrations comparable to those in eutrophic estuaries globally [38], with a 40-and 800-fold increase in medium and high treatments, respectively. The large gradient of environmental conditions encompassed within this field study ensured our estimates of MPB productivity fell within ranges reported within intertidal habitats both in New Zealand and worldwide. For example, NPP estimates have previously been described between 500 and 3000 µmol O 2 m −2 h −1 within emerged habitats [25,47] and between −3000 and 4000 when measured continuously over several tidal cycles [48].
MPB productivity measurements were spatially and temporally variable and, even after normalising for background site values (i.e., NPP CN and GPP CN ), both emerged and submerged MPB production were independent of nutrient enrichment (Table A5; Figure 3). This was unexpected considering each site was considered to be nitrogen-limited (within control plots) according to the calculated N:P porewater ratios ( Figure A2) [45,46]. The negligible role of sediment nutrient enrichment may be explained by the dilution of plot effects from the suspension and removal of MPB with bedload transport [49], and/or the potential oversupply of nitrogen relative to phosphorus, switching the system from being predominately nitrogen to phosphorus limited. However, the absence of a response from MPB to stoichiometrically balanced nutrient additions have previously been observed [30]. This may occur if nutrients were not limiting, for example, as a product of both water-column assimilation and sediment porewater nutrients, the replenishment of depleted nutrients through advective flushing (by macrofauna and physical processes [50,51]), and/or through intense grazing pressure reducing MPB biomass and therefore reducing N demand [31]. Additionally, light limitation can prevent MPB growth [32,52], which while plausible for sites in this study with high water column turbidity (MNK, RAG), is unlikely to prevent a response by MPB at the least turbid sites. We therefore postulate that a combination of factors may have prevented an observed increase in MPB growth in this study, such as high bedload transport, carbon supply, light availability, phosphorus limitation, or grazing activity and therefore the enrichment of sediment nitrogen had no observable effect on MPB production.
In contrast to nutrient enrichment effects, site turbidity had a substantial influence on primary production estimates. For example, NPP was negative at the most turbid sites and increased with decreasing site turbidity. Although comparisons of GPP across sites did not show as strong a relationship with water column turbidity, submerged GPP was greatest at two of the least turbid sites (TAU-O, TAU-T). These results are further supported by the DistLM analysis where 40 and 17% of submerged NPP and GPP variability, respectively, could be explained by site turbidity (Table 6). While there was still a large proportion of unexplained variation, this is most likely an artefact of dynamic variations in sediment biogeochemistry caused by microbial and macrofaunal activities and detrital inputs [43,53,54] which will all contribute to small scale variations between plots. Nevertheless, the strong relationship with site turbidity appears to dominate benthic primary production estimates over time by occluding the light required by MPB for photosynthesis. These results highlight the effectiveness of incorporating a natural gradient of water column turbidity as a way of understanding the complexities of multiple stressors on ecosystem functioning in real-world settings.
As the tide recedes and exposes the sediment surface, MPB production becomes unconstrained by light limitation, confirmed by our results showing decreased effects of light on MPB production during emerged tidal periods. For example, PAR Site explained just 11 and 3% of emerged NPP and GPP variability, respectively. As the effects of light decreased, the influence of sediment characteristics and, to a lesser extent, temperature, increased (GPP, mud 14%; NPP, temperature 4%; Table 6). However, considering the large influence temperature can have on MPB metabolism and previously reported production estimates [23,52,55], we postulate the amount of variability explained by temperature may increase if sampled over a larger temperature gradient (13.6-18.7 C during emerged periods in this study). Sediment characteristics, in contrast, can modify primary production estimates through several direct and indirect pathways. These include alterations in the penetration depth of light into the sediment, changes in solute transport in less permeable sediments, and a reduction in the metabolic status of MPB [56]. Moreover, sediment properties can also alter MPB community composition, resulting in different assemblages between sites [57,58]. While this study incorporated the effects of these environmental differences (including turbidity responses) and the effect of long-term nutrient enrichment on MPB production, characterisation of MPB assemblage responses presents an interesting avenue for future research.
The disparity in environmental conditions experienced during emerged vs. submerged periods is likely to be a contributing factor to the observed difference in photosynthetic efficiency between tidal states. During periods of submergence, few differences in photosynthetic efficiency were observed between sites, suggesting potential photoadaptation to lower light availability. Once released from light limitation during emerged periods, photosynthetic efficiency was significantly lower when compared to submerged tidal periods, and this was consistent through time ( Figure 5). During periods of emergence, MPB efficiency can be reduced through high light and UV-B exposure, temperature extremes and desiccation stress (and the subsequent changes to salinity) [55,59,60]. However, a consistent difference in photosynthetic efficiency was not observed between the summer and winter sampling periods and no change in sediment porosity (and therefore desiccation) was detected over each emerged period (Table A2), suggesting these factors were unlikely to significantly reduce photosynthetic efficiency in this study. Alternatively, protective mechanisms such as diel migration are likely to have altered photosynthetic efficiency during emergence [61].
At the three most turbid sites, the contribution of emerged productivity was significantly greater than submerged (Tables 4 and 5), supporting a previous study suggesting the importance of emerged MPB production increases with site turbidity [25]. In addition, the significant contribution of emerged primary production to total productivity is highlighted by the dominance of net heterotrophy at highly turbid sites during submergence, compared to predominately net autotrophy at the least turbid sites ( Figure 6). This suggests highly turbid sites may rely on emerged periods (which were net autotrophic across all sites) to sustain and support benthic primary production. This is in agreement with studies from around the world reporting the restriction of benthic primary production to emerged periods only, owing to light limitation through high turbidity [21,22,62]. Despite the maintenance of MPB productivity over emerged periods, the higher production during emergence does not always fully compensate for the lower production during submergence, as evidenced by the consistently higher MPB production estimates during submerged compared to emerged tidal periods at the least turbid sites (Figures 4 and 5). Increased heterotrophy has further consequences for coastal ecosystem functioning via the associated release of inorganic nutrients through increased respiration (rather than assimilation to support productivity during periods of net autotrophy). This ultimately can have substantial implications for coastal ecosystems through indirectly facilitating the removal of nutrients in oligotrophic systems or contributing to accelerating water column primary production in eutrophic systems. Consequently, reduced MPB productivity during submergence not only has direct implications on the supply of labile carbon and therefore coastal food webs [6] but can cascade onto reductions in the capacity to moderate pollutants, changes in trophic structure and alterations in nutrient cycling and transformation pathways [10,63,64]. In combination, these changes feedback within the system to modify ecological interaction networks and push ecosystems closer towards tipping points [65][66][67].
Through coupling measurements of submerged and emerged primary production, this study highlights the multifaceted nature of intertidal habitats where productivity is moderated by a complex network of ecosystem interactions. However, it is clear that reductions in seafloor light climate leading to any decline in primary productivity by MPB will have a cascading influence on the health and functioning of coastal habitats. Our study highlights the importance of low tide MPB production in turbid estuaries and therefore as a contributor to the system's resilience to elevated turbidity. Intertidal areas are already vulnerable and it is estimated that there has been a 16% loss of intertidal habitats globally [68], with this expected to increase in the future [69]. As emerged tidal habitat becomes restricted and lost, the ecological resilience that these areas were providing is also lost, having significant implications for the effective management of these valuable ecosystems.        Table A5. A repeated measures PERMANOVA using Euclidean distance was used to detect treatment effects on submerged and emerged primary production estimates. The analysis had treatment (3 levels), site (6 levels) and time (2 levels; sub T1 and T3; em T2 and T3) as fixed factors and replicate plot nested within treatment. Significant effects (p < 0.05) are given in bold. Main effects were not considered if interaction effects were significant and post hoc pairwise tests were undertaken to identify where differences occurred.