Biodiversity breakpoints along stress gradients in estuaries and associated shifts in ecosystem interactions

Denitrification in coastal sediments can provide resilience to eutrophication in estuarine ecosystems, but this key ecosystem function is impacted directly and indirectly by increasing stressors. The erosion and loading of fine sediments from land, resulting in sedimentation and elevated sediment muddiness, presents a significant threat to coastal ecosystems worldwide. Impacts on biodiversity with increasing sediment mud content are relatively well understood, but corresponding impacts on denitrification are uncharacterised. Soft sediment ecosystems have a network of interrelated biotic and abiotic ecosystem components that contribute to microbial nitrogen cycling, but these components (especially biodiversity measures) and their relationships with ecosystem functions are sensitive to stress. With a large dataset spanning broad environmental gradients this study uses interaction network analysis to present a mechanistic view of the ecological interactions that contribute to microbial nitrogen cycling, showing significant changes above and below a stressor (mud) threshold. Our models demonstrate that positive biodiversity effects become more critical with a higher level of sedimentation stress, and show that effective ecosystem management for resilience requires different action under different scenarios.


Results
Quantile regressions (90%) of biodiversity measures and sediment mud content showed positive relationships at very low levels of mud (~4% on average), after which the relationships became negative (Fig. 1). DEA increased with both mud and organic content, and although mud and organic content are strongly correlated (Pearson's R = 0.69, p < 0.001), the sandy (<4% mud) and muddy (>10% mud) data subsets showed overlapping values of organic content as well as other model response variables (Figs. 2 and 3). The conceptual basis model used for our structured equation model analyses demonstrates the high potential for multiple direct and indirect controls on DEA, and the highly connected system that contributes to variability in DEA (Fig. 4). Interaction networks were significantly different in sandy compared with muddy systems (Fig. 5), and all response variable R 2 values decreased in the muddy model (Fig. 5). Most notable was the reduction in significant positive network linkages (from 9 to 6) and increase in the number of significant negative network linkages (from 1 to 3) from the sandy to muddy model (Fig. 5). Sediment organic content changed from having   Box plots for ecosystem components (a) mud content, (b) organic content, (c) Chl a, (d) number of taxa, (e) macrofaunal abundance, and (f) number of large bivalves in the sandy (<4% mud, n = 75) and muddy (>10% mud, n = 57) datasets. Boxes show 25 th and 75 th percentiles, whiskers the 10 th and 90 th percentiles and black dots the 5 th and 95 th percentiles. Solid line is median and dashed line is mean.

Discussion
Our dataset was consistent with other studies demonstrating changes or breakpoints in biodiversity measures with increasing mud content 18,26,29,38 . The dataset showed a high degree of variability as expected in real-world ecological systems [47][48][49] ; for example, macrofauna-sediment relationships showed factor-ceiling-like relationships (Fig. 1). This high degree of variability can make assessing or detecting changes in important ecological relationships and interactions difficult, however SEMs and interaction networks help us disentangle the relationships present among this variability and understand the underlying ecosystem processes 50 . The differences in ecological interaction networks observed for <4% and >10% mud content demonstrated how potentially important ecological relationships can change with the muddying of intertidal estuarine habitats. Our results are consistent with other work showing ecosystem changes associated with increases in sediment mud content including; a threshold response at 20% mud 28 , reductions in macrofaunal traits based indices above 10% mud 51 , and functional shifts in the relationship between a key consumer and the abundance of microphytobenthos above and below 15% mud content 49 .
Although muddy sediments are typically organic rich, sediment mud content may affect ecosystem functioning, especially denitrification, in a concentration-dependent or threshold type manner due to interactions between mud and OM and their combined effects on macrofauna. Strong correlations among these drivers of DEA (i.e. mud and organic content) make it hard to decipher whether effects of changing sedimentary environment are due to physical effects caused by higher proportions of finer sediments, or chemical effects such as supply of bioavailable nitrogen and carbon from breakdown of organic matter. Muddy sediments tend to have high organic content and therefore pore water nitrogen concentrations (from organic matter breakdown), but slower rates of pore water advection compared with sandy sediments 52 . For microbial nitrogen cycling this means that although there is higher supply of inorganic nitrogen in muddy sediments, rates of nitrogen transformation are limited by the movement of solutes throughout the sediment, whereas in sandy sediments there is faster pore water flow and solute transport, but rates of nitrogen processing are limited by lower solute concentrations 12 . Our models reflect this; the influence of mud versus organic content on DEA differed in our sandy (<4% mud) and muddy (>10% mud) datasets. In sandy systems, organic content is an important driver of DEA (and other ecosystem response variables), suggesting that DEA is more limited by the quantity of inorganic nitrogen substrates. Whereas our muddy sediment SEM suggests that DEA is more limited by the supply of resources to specific sediment horizons where they can be transformed (due to the slow speed of diffusion through the smaller pore spaces in cohesive muddy sediments). The influence of sediment organic content on DEA is less in the muddy model, and the effect of macrofaunal abundance on DEA remains an important linkage, despite the lower path coefficient and much lower explained variance of macrofaunal abundance as a response variable (R 2 = 0.21). It follows that ecosystem components that enhance sediment mixing (i.e. bioturbation), such as macrofaunal abundance, may become relatively more important in muddy sediments because the bioadvection of solutes facilitates the supply of resources to DEA (and other biogeochemical processes).
Bioturbation is an important ecosystem function carried out by benthic macrofauna that is known to significantly enhance microbial nitrogen processing 17 . The effect of bioturbation on biogeochemical fluxes is enhanced by the abundance of key species (large taxa, such as A. stutchburyi and M. liliana in the present study, that mix and www.nature.com/scientificreports www.nature.com/scientificreports/ mobilise sediments via their movements and feeding behaviours), community abundance, species richness, and functional diversity 5,53-56 . These models suggest that the link between macrofaunal abundance and DEA represents a critical positive relationship (and link between other important influences on DEA including macrofaunal diversity and microphytobenthic biomass (Chl a)), and that the strength of this link is reduced with increasing stress (mud content) and the associated changes in the overall ecosystem interaction network.
These network models allow us to characterise how ecosystems respond to increasing stressors at a system level by helping to identify the mechanisms behind ecosystem function response to perturbation, as well as understand the consequences of how further anthropogenic pressure may alter the system, its interactions and ecosystem functions 57 . Our interaction network analyses show how changes in an ecosystem, such as shifts in macrofaunal community, may have effects on other parts of an ecosystem, or ecosystem functions, and that shifts may have differing effects under different scenarios or stress regimes. For example, overharvesting of shellfish (large bivalves) may have several indirect effects on microbial nitrogen cycling (DEA), but these effects are likely to be more severe in ecosystems with elevated mud content. Similarly, loss of biodiversity and abundance of macrofauna, will have differing effects on ecosystem dynamics, and therefore microbial nitrogen processing, in estuaries that have high sedimentation levels (high mud content) than those that are sandier. This result may www.nature.com/scientificreports www.nature.com/scientificreports/ be indicative of how important ecological relationships may change with increases in other stressors, however further analysis including measures of other stressors (e.g. nutrient enrichment) and ecosystem functions (e.g. primary productivity) is required to test this. Comparing the two models reveals that higher mud content resulted in a reduction in the amount of explained variability in ecosystem components (including DEA), and this can be symptom of increased ecosystem stress [58][59][60][61] . Increased variability in ecosystems can also be indicative of an approaching tipping point or regime shift [62][63][64] , thus our approach could be used as an early warning indicator as well as a tool for mitigating such changes.
It should be noted that many healthy estuaries are likely to fall within the 4-10% mud content range not analysed in our study, and <10% mud content would ordinarily not be considered impacted by sedimentation. We do not have a sufficient number of data points within the 4-10% mud content range to build a SEM, but we would anticipate some differences in linkages compared with the low mud model. Sediment stability (cohesive to non-cohesive), biogeochemical processes and macrofaunal communities and their inter-relationships change substantially across this range of mud content 25,29,65 and we would expect the interaction network to reflect this accordingly. With the data we do have between 4-10% mud content Pearson's correlation coefficients revealed strong positive correlations between large bivalves and DEA, and Chlorophyll a and DEA, suggesting that biological and macrofaunal ecosystem components are again key drivers of functionality. Without a larger dataset we cannot speculate on whether key functional linkages would be lost within this range but it is worthy of further study.
Investigating ecosystem resilience and thresholds is increasingly important in an era of unprecedented anthropogenic change. Detecting and responding to major ecosystem changes before they occur is a critical challenge for ecologists and managers of coastal ecosystems 64,66,67 . However, the way to quantitatively assess resilience and detect thresholds (from a management perspective) is less clear and requires more research 68 . Our investigation has demonstrated that at different levels of a stressor (sediment mud content) the architecture of ecosystem interaction networks changes with consequences for ecosystem function. Our structural equation model analyses support previous experimental work demonstrating the interdependence of estuary microbial nitrogen processing on both physical (sedimentary environment) and biological (macrofaunal community and biodiversity) ecosystem variables 5,12,13 . Furthermore, we show evidence for loss of resilience of microbial nitrogen cycling with increasing sedimentation stress through changes in interaction networks, increased negative effects on potentially important network linkages (especially biodiversity related), and increased dependence on these positive biodiversity effects. Our models showed that the mechanisms of ecosystem components acting on DEA are different under different levels of sedimentation stress, and therefore managing for ecosystem resilience to nutrient enrichment requires different action in different scenarios.
This study confirms that there is high interdependence of microbial nitrogen processing on other ecosystem components and their interactions, meaning that resilience to nutrient enrichment will be highly responsive to ecosystem changes such as increasing sediment mud content. Interaction networks reveal that DEA is more tightly linked, or dependent on, the macrofaunal community in muddy sediments than in sandy sediments. In agreeance with many other studies, this means that the effects of ecological relationships, especially those involving biodiversity components, are context dependent, and more imperative to ecosystem resilience when levels of a stressor are higher [69][70][71] . This research is an example of using data from site-specific studies to scale up and increase the generality and broadness of understanding of how ecological relationships and interaction networks operate over greater environmental gradients and changes in ecosystem components. Preservation of benthic macrofaunal community biodiversity is essential to mitigate the effects of increasing stressors on nitrogen cycling and ecosystem resilience, and this has greater importance in ecosystems already under sedimentation stress. Understanding these general patterns and how they change across environmental gradients, and/or with increasing stressors provides a way forward for management and conservation when working with ecological datasets that are inherently heterogeneous.

Methods
Data were compiled from two published (the un-manipulated/ambient plots in Douglas, et al. 69,72 ) and three unpublished surveys 73,74 carried out in five northern New Zealand estuaries with extensive intertidal areas, and catchments with variable land use (Fig. 6, Supplementary Table S1). In total, the dataset consisted of samples collected from 171 plots (1 m 2 ) during austral summer months (December-March) between 2013 and 2018. The estuaries contained overlapping gradients in abiotic (especially sediment properties) and biotic ecosystem components allowing data to be pooled and analysed for patterns across estuaries.

Denitrification enzyme activity-the primary response variable of interest. Methods for
Denitrification Enzyme Activity (DEA) sample collection and assays are described in full in Douglas, et al. 69 , but briefly, 5 replicate sediment cores (0-5 cm depth, 5.3 cm dia.) were taken randomly from each plot and pooled. Unfiltered seawater was collected from each site, and water and sediment were stored on ice for transport to the laboratory where they were stored at 4 °C. Assays were conducted at room temperature (20 °C) within 48 h of sample collection using the chloramphenicol amended acetylene inhibition technique 41,75 . Assay slurries contained 30 mg L −1 C as glucose and 10 mg L −1 N as KNO 3 . Headspace gas samples (6 mL) were collected from each assay over a 2 h time course (10,30,60,120  www.nature.com/scientificreports www.nature.com/scientificreports/ macrofaunal invertebrates >500 µm longest dimension were sampled using 13 cm dia., 15 cm depth sediment cores (n = 2 per plot in Kaipara Harbour, n = 1 per plot at all other sites) 76 . For measurements of sediment mud content (Mud), organic matter was first removed from samples using 10% hydrogen peroxide, then samples were analysed using a Malvern Mastersizer 2000 77 to derive the proportion of mud (particle sizes <63 µm). Sediment samples were dried to a constant mass (60 °C, 48 h), then organic content (OC) was determined by weight loss on ignition (550 °C for 4 h) 78 . Microphytobenthic biomass (Chlorophyll a (Chla)) was determined after extraction of pigments from freeze-dried sediments using 90% acetone, then a Turner Designs 10-AU flourometer was used to measure fluorescence 79 . Preserved (50% isopropyl alcohol) macrofauna were stained with Rose Bengal and in the laboratory, all organisms were counted and identified to the lowest possible taxonomic level (usually species). Macrofaunal community variables (biodiversity measures) included total abundance (N), number of taxa (S), and number of the large bivalve species Macomona liliana and Austrovenus stutchburyi (LB; considered separately because of their substantial influence on ecosystem function on New Zealand sandflats 26,29,80 ).

Low versus high mud datasets.
Our key question was how ecological interaction networks change with increased sediment muddiness and how this may affect DEA. This required the development of SEMs for sandy (low mud) and muddy situations. We created two datasets (sandy and muddy), given the following constraints. SEMs require a large number of data points, with the minimum number of data points determined by the number of nodes present in an initial conceptual basis model (see below). Any split of the data into sandy vs muddy required at least 55 data points in each subset. Further, we endeavoured to create two data sets with an approximately equal number of data points in each, as this was the most objective way to compare the significance of paths in sandy and muddy SEMs.
Bivariate plots of biodiversity measures with sediment mud content as the predictor variable were examined to investigate thresholds in mud content that could be indicative of changes in ecosystem interaction networks (Fig. 1). This was done using quantile regression and breakpoint analyses using packages 'quantreg' 81 and 'segmented' 82,83 in R 3.5.1 84 . Breakpoints in sediment mud content (the predictor variable) for each biodiversity measure (number of species, number of individuals, and number of large bivalves) were determined using breakpoint analyses of the 90 th percentile (τ = 0.9) quantile regression models (Fig. 2). The average breakpoint was determined to be 4% mud, which was then used to create a low mud content dataset (i.e., "sandy", <4% mud, n = 75). A higher mud content dataset was created with values >10% mud content (i.e., "muddy", n = 57). This value is based on previous work showing that DEA responds negatively to stress when there is >10% mud content (but not <10%) 72 , indicating potential changes in ecosystem interaction networks associated with DEA may occur around this point. Having disjunct datasets (i.e., <4% mud for sandy, and >10% mud for muddy) allowed us to present coherent results for mud levels above and below the 4-10% mud range, which literature suggests may encompass other ecologically significant thresholds 28,51,72,85 . Specifying a precise breakpoint position within this range was not required. Given our analytical goals and constraints (e.g., minimum dataset size, relatively balanced datasets), our approach was appropriate. Each dataset contained a minimum of four samples from each estuary, except for Mahurangi where there were only samples with >10% mud content. This minimised any confounding effects of 'estuary' on results.