Lagoon Biogeochemical Processing is Reflected in Spatial Patterns of Sediment Stable Isotopic Ratios

The spatial analysis of biota, particulate organic matter, and sediments for stable isotopes of carbon (δ13C), nitrogen (δ15N), and sulfur (δ34S) have proved useful for identifying patterns in productivity, nutrient pollution, and relationships between biological and physiochemical variables at the local and global scales. Yet such approaches are rarely applied to studies of lagoon or estuarine metabolism. Focusing on Bahía San Quintín, a heterotrophic seagrass-dominated lagoon on the Pacific coast of Baja California, México, we report on spatial patterns in surficial sediment CNS stable isotopic ratios as tracers of lagoon biogeochemical function. Stable nitrogen isotopes highlighted potential spatial variability in the balance between denitrification and nitrogen-fixation within the lagoon and identified an association between elevated δ15N levels and oyster culture, suggesting that oyster presence may be enhancing N2 production. Spatial patterns in δ34S covaried with sediment particle size, underlining the importance of sediment texture in determining the depth of sub-oxic-anoxic redox zones. Sediment carbon stable isotope ratios highlighted the lack of incorporation of seagrass carbon into seagrass meadow sediments, thus emphasizing the importance of phytoplankton or microphytobenthos for carbon accumulation in seagrass meadows. This report highlights the value of sediment isotopic values in corroborating spatial patterns in estuarine metabolism or macronutrient processing identified from chamber or flux-based studies. Stable isotope mapping can provide a useful addition to assessment of estuarine metabolism, or act as a stand-alone tool for generating hypotheses, identifying the influence of spatial gradients, and/or suggesting prime locations for investigation of microbial abundance or function.


Introduction
In recent years, the analysis of biota, particulate organic matter, and sediments for stable isotopes of carbon (C), nitrogen (N), and sulfur (S) has emerged as a common and powerful approach for assessing estuarine, nearshore, and open-ocean ecosystems [1,2]. In nearshore environments, δ 15 N variation has been linked to marine vs. terrestrial N sources, upwelling, N-fixation, denitrification, and wastewater pollution [3,4]. Variation in δ 34 S has been shown to reflect salinity, sulfate reduction, and benthic organic matter sources, and δ 13 C variation has been linked to rates and patterns of 2 of 19 primary production, depth, and aqueous CO 2 signatures [3,5]. However, despite the demonstrated ability of landscape and regional-scale stable isotope analysis to uniquely reveal spatial gradients in ocean productivity [6], nutrient pollution plumes [7], and nutrient recycling [8], studies have only rarely applied this approach to coastal estuaries and lagoons [9]. Here, we focus on the ability of spatial patterns in CNS sediment isotope ratios (δ 15 N, δ 34 S, δ 13 C) to provide insight into patterns of macronutrient processing across a lagoonal landscape, including potential effects of seagrass presence and oyster aquaculture on biogeochemical cycling.
Rates of N cycling are often elevated in coastal marine habitats due to high organic matter inputs coupled with anaerobic respiration processes that use nitrate as a terminal electron acceptor [10]. Until recently, it was believed that denitrification, or the bacterially-mediated reduction of nitrate to N 2 gas, dominated coastal N-processing [11]. Denitrification, which can return anthropogenically-derived reactive N to non-reactive N 2 gas and therefore buffer the open ocean from the impacts of human-derived N pollution [12], has been particularly well-studied in the context of its role as an ecosystem service of societal value [13]. In addition, denitrification in the three ocean pelagic oxygen minimum zones is thought to be responsible for a large percentage of the total oceanic loss of N, and has therefore been a focus of investigation into global marine biogeochemical cycles [14]. However, with the widespread usage of isotopic tracer techniques [15], other pathways of N cycling have emerged as common constituents of the oceanic N cycle, including N-fixation [16]. Although N−fixation, or the microbially-mediated conversion of N 2 gas to reactive N (NH 4 + ), has traditionally been regarded as a minor coastal N source [17,18], recent studies have highlighted its importance [19,20]. N-fixation has been found particularly in oligotrophic coastal marine environments [21], tropical and sub-tropical waters inhabited by diazotrophic photoautotrophic cyanobacteria [22], and in association with seagrasses, whose exudates provide a labile C source to support N-fixation by heterotrophic bacteria [23]. Although the balance of denitrification vs. N-fixation in most coastal lagoons remains unresolved, stable N isotopes can be useful for discriminating N-fixation (which introduces bioavailable N at isotopic values ≈ 0% ), and denitrification which consumes lighter N at a faster rate, leaving the NO 3 pool substantially enriched in 15 N [5]. Sulfate is the most important electron acceptor for the anaerobic oxidation of organic matter in coastal sediments [24], due to the abundant pools of sulfate in seawater compared with other potential terminal electron acceptors [25]. Sulfate reduction is thus a key process in sub-oxic and anoxic marine environments. Sulfides are produced from sulfates and tend to accumulate dissolved in porewater, as gaseous sulfide, or precipitate with iron to form iron sulfide minerals (e.g., pyrite). The global average value for δ 34 S in seawater is +20.7% [26], but sulfate is fractionated during cycles of sulfate reduction and disproportionation, resulting in sulfide increasingly depleted in 34 S [27]. Consequently, sedimentary sulfide has a δ 34 S signature of −15 to −25% , which is preserved without fractionation when assimilated by plants [28]. Sulfur stable isotopes thus act as a tracer for anaerobic decomposition, in addition to benthic organic matter sources in coastal food webs.
Stable C isotopes in marine sediments or particulate organic matter often act as productivity tracers because of the partial CO 2 limitation that occurs in highly productive aquatic environments. Photosynthesis discriminates against heavier C isotopes, but this photosynthetic fractionation is lower during periods of high growth rates and low aqueous CO 2 concentrations [29]. However, in highly productive nearshore environments, including in seagrass meadows, light availability often acts as a primary control on production such that C stable isotopes vary along depth gradients [3,30]. In addition to depth and productivity, δ 13 C in estuarine sediments can also reflect organic matter sources (e.g., seagrass, phytoplankton, algal, terrestrial) due to divergence in photosynthetic pathways that are associated with different discrimination factors and C sources (CO 2 vs. HCO 3 − ) [31].
This combination of naturally abundant stable isotopic tracers is especially relevant to understanding patterns of macronutrient cycling in seagrass-dominated coastal lagoons. Seagrasses rank among the most productive marine primary producers, and as ecosystem engineers that modify currents, trap sediments, and constrain benthic microbial communities, seagrass meadows modulate the mineralization and regeneration of organic matter [32]. Although net effects differ by species, seagrasses are known to modify biogeochemical conditions by exuding labile C which stimulates microbial processes and sulfate reduction, in addition to exuding oxygen, which promotes root zone oxidation [33]. In highly productive seagrass meadows, diel cycles in water column dissolved inorganic carbon (DIC), O 2 , and pH may occur, and seagrass beds may also act as both a source and a sink for C depending on the bed age and season [32]. Seagrass C and S stable isotopic ratios tend to be distinct from other marine primary producers because they readily assimilate bicarbonate [34] and experience sulfide intrusion [35], which result in heavier C and lighter S stable isotopic compositions, respectively. Although previous work has used this distinct isotopic profile to trace seagrass-derived organic matter in coastal sediments and food webs [36] and assess seagrass health relative to nutrient pollution and sulfide intrusion [35,37], few studies have focused on understanding how seagrass presence modifies patterns of sulfate reduction, N fixation, or denitrification at the landscape scale [38].
Another factor thought to shape coastal nutrient cycling is the establishment of shellfish aquaculture. As an important and growing industry [39], cultured shellfish may contribute to altered macronutrient cycles in coastal lagoons through depletion of phytoplankton, creation of biogeochemical hotspots in association with bio-wastes, and nutrient bioextraction associated with shellfish harvests and denitrification [13,40]. Although shellfish aquaculture of non-native species is often considered a disturbance from the perspective of sustainable management of Pacific coast estuaries [41], the reintroduction of native oysters and mussels is thought to support improved water quality in Atlantic estuaries [42]. Oyster restoration has also been associated with other beneficial services such as improved water clarity, habitat provision for fisheries, and wave attenuation [43,44]. Although previous studies of long-line shellfish aquaculture at low stocking densities have generally noted only minor habitat impacts [45], effects of aquaculture on lagoon nutrient cycling may extend further than culture racks. Therefore, impacts of shellfish aquaculture may be assessed through estuary-wide, spatially distributed sampling.
Here, we report on spatial patterns in surficial sediment stable isotope ratios (δ 15 N, δ 34 S, δ 13 C) in Bahía San Quintín (BSQ), an inverse lagoon on the Pacific Coast of Mexico, as tracers of lagoon biogeochemical function. The absence of terrestrial inputs due to highly episodic riverflows (once per decade) and intense groundwater overdraft that has reversed the flow of shallow groundwater [46] makes BSQ an ideal location to study oceanic-driven nutrient processing. Stable N isotopes were used to identify sedimentary signatures indicative of oceanic N (δ 15 N = 8-10; [47]), vs. contributions from denitrification or N-fixation, which would be reflected as enriched or depleted δ 15 N signatures, respectively [48]. S (δ 34 S) stable isotopes were used as tracers of sulfate reduction, which occurs in concert with anaerobic sediment respiration. Because BSQ is net heterotrophic [49], and anoxic respiration has been estimated to account for slightly less than half of the respiration of BSQ [50], sulfate reduction is likely an important tracer of benthic metabolism. Finally, C isotopes (δ 13 C) were evaluated relative to organic matter sources, aqueous CO 2 sources, and spatial patterns of productivity and respiration. We used these stable isotopic tracers to focus on the following objectives: (1) to identify potential spatial patterns in N-fixation, denitrification, and S reduction, and relate these patterns to known spatial gradients in productivity, depth, and benthic sediment characteristics, and (2) to identify how seagrass cover and oyster culture are associated with patterns of altered biogeochemical function, with the ultimate goal of advancing the use of sediment isoscape data to inform studies of lagoon biogeochemical dynamics.

Materials and Methods
Research was conducted at BSQ ( Figure 1A), a hypersaline lagoon located on the northwestern coast of Baja California (30.4152 • N, −115.9541 • W). The lagoon has an area of 42 km 2 , is covered by extensive seagrass beds, and experiences semi-diurnal mixed tides that range from 1.0 m (neap) to 2.5 m (spring). The structure of the lagoon is defined by the San Quintín volcanic field, which includes eight volcanic complexes and two isolated volcanos (Mt. Mazo and Isla San Martín), which date to Pleistocene and Holocene in age [51]. The lava fields formed as subaqueous volcanos which emerged as islands, and then acted as traps for sand transported by longshore drift [52]. The western boundary of the lagoon is a large tombolo separating two volcanic cones [53], while the sand spit that forms the southern limit of the lagoon has been built by wave refraction and northerly drift in the lee of the outer spit and cones [54]. southern limit of the lagoon has been built by wave refraction and northerly drift in the lee of the outer spit and cones [54]. Rainfall at BSQ averages 143 mm y −1 ; yearly evaporation averages 1436 mm y −1 [55]; and the lagoon experiences surface water inflow on average every 10 y, and groundwater inflow on average every 4 y [46]. Thus, the lagoon is persistently hypersaline throughout the year, with salinity levels increasing from the oceanic inlet towards the heads [49]. Tidal exchange occurs through a single mouth that divides into two arms: the eastern arm is known as Brazo de San Quintín, and the western arm is known as Bahía Falsa (BF). The depth of the lagoon is less than 2 m, except in channels where depths range from 5 to 8 m [56]. Cultivation of the pacific oyster, Crassostera gigas, was introduced to Rainfall at BSQ averages 143 mm y −1 ; yearly evaporation averages 1436 mm y −1 [55]; and the lagoon experiences surface water inflow on average every 10 y, and groundwater inflow on average every 4 y [46]. Thus, the lagoon is persistently hypersaline throughout the year, with salinity levels increasing from the oceanic inlet towards the heads [49]. Tidal exchange occurs through a single mouth that divides into two arms: the eastern arm is known as Brazo de San Quintín, and the western arm is known as Bahía Falsa (BF). The depth of the lagoon is less than 2 m, except in channels where depths range from 5 to 8 m [56]. Cultivation of the pacific oyster, Crassostera gigas, was introduced to BF in the 1970s, because BF has higher phytoplankton concentrations [49], and up to 5000 tons are produced there annually [57]. Coastal upwelling associated with the California Current System is considered the main source of nutrients to BSQ, and upwelling tends to be most intense in the spring-summer period when northwesterly winds are strongest.
Surface sediment grabs were collected from BSQ in November 2016 using a Van Veen sediment grab sampler (69 samples), and in July 2018 using a petite Ponar grab sampler (15 samples), and processed for bulk stable isotopic ratios, organic and inorganic C content, N content, particle size distribution, loss on ignition, and bulk density. Benthic cover of seagrass was recorded using a submerged camera (FDR X3000, Sony, Tokyo, JP) or videography system (SeaDrop, 550 Seaviewer, Tampa, FLA). Distribution of seagrass was digitized using maximum likelihood supervised classification of 4-band Planet imagery (2 July 2017) taken at low tide using ENVI 5.4 (Harris Geospatial, Broomfield, CO, USA, 2016), and distribution of oyster racks was manually digitized using low-tide NIR drone imagery (eBee senseFly, 0.07 m resolution, collected 29 November 2016). A cost-distance function was used to estimate the product of the least accumulative cost-distance function calculated as a raster in ArcGIS 10.5, using the product of the distance to the mouth and friction categories based on depth.
Due to instrument availability, stable isotope data was provided by three different isotope ratio mass spectrometers (IRMSs), including a CHNOS Elemental Analyzer interfaced to an IsoPrime100 IRMS at UC Berkeley, an Elementar Pyrocube interfaced with an Elementar IsoPrime100 IRMS at the Academy of Natural Sciences of Drexel University, and a COSTECH Elemental Analyzer coupled to a Thermo DeltaV Advantage IRMS at the Laboratorio de Isótopos Estables (LIE) at CICESE. Stable isotope values are expressed relative to the isotopic ratio of a standard following Formula (1): where a R X represents the ratio of the less common to the more common isotope in the sample and a R Std represents the ratio of the less common to the more common isotope in the standard reference material (air, Vienna Pee Dee Belemnite, and Vienna Canyon Diablo Troilite for δ 15 N, δ 13 C, and δ 34 S, respectively). Mean sample differences for replicates were 0.8% for δ 13 C, 0.84% for δ 15 N, and 0.48% for δ 34 S. For samples run on different instruments, the mean sample difference for replicates was 1.09% . Samples analyzed for C isotopes were acidified to remove inorganic C and samples analyzed for N isotopes were not pre-treated. Samples processed for S isotopes were triple-rinsed using deionized water [58] to remove sulfate, which could have removed oxidized iron monosulfides from the sample. Thus, the reported δ 34 S values reflect sulfur of organic matter and pyrite only. Potential sources of organic matter to the lagoon sediments, including Zostera marina, Ulva spp., Crassostrea gigas feces and pseudofeces, and marsh macrophytes, were also collected and processed for stable isotopes (Table 1). Particle size distribution was measured on sediment samples pre-treated with 30% H 2 O 2 to remove organic material [59] then introduced into a Beckman Coulter LS-13-320 system recording particle size distributions from 0.04-2000 µm with 117 bins, and run in triplicate. Average particle size distributions were post-processed with Gradistat.v8 software [60], including bin aggregation to texture classes, and statistical description. Sediment textures were defined per [61]. Particle size (diameter, d) values were converted to phi (ϕ) units using Formula (2): The full logarithmic method of moments was used to calculate particle size mean, sorting, skewness, and kurtosis. Sediment samples were processed for bulk density and loss on ignition by drying a known volume of sediment to constant weight and combusting for 4 h at a temperature of 550 • C [62].
Sediment stable isotope ratios and associated sediment data (e.g., elemental ratios, grain size distribution) were visualized through creation of contour maps using ordinary kriging, utilizing a circular semivariogram model and assuming isotropic conditions. For data analysis, depth and residence time were estimated for each sample point using previous hydrographic surveys and model outputs [63,64]. Statistically significant hot and cold spots in the stable isotope data were identified using Anselin Local Moran's I statistic, and differences in isotopic composition were identified between the two arms of the lagoon using the Mann-Whitney-Wilcoxon test. Regression analysis (Mann-Kendall) was used to identify significant spatial trends in sediment isotope data based on latitude and the product of the least accumulative cost-distance function. Isotopic signatures were compared graphically and statistically for magnitude and variance (using Mann-Whitney-Wilcoxon and Fliger-Killeen tests, respectively) for: (1) sediment samples removed from seagrass beds vs. from bare sediments, (2) sediments collected from within or adjacent to seagrass meadows (<100 m) vs. locations distant from sea grass meadows, and (3) sediment samples taken within (<100 m) areas of oyster culture vs. more distant. Stable isotope ratios were modeled using a backward stepwise approach as a function of sediment characteristics and grain size metrics using principal least squares (PLS) regression to account for the high collinearity among predictors. Collinearity of predictors was visualized through creation of a non-parametric co-variance matrix, ordered using the first principal component. Geospatial analyses were completed using ArcGIS ver. 10.5, Redlands, CA, and statistical tests were conducted in R 3.6.1 [65] using packages pls [66], plsVarSel [67], Hmisc [68], corrplot [69], and Kendall [70].

Results
Here we report bulk stable isotope ratios (N, S, C) of surface sediments to identify spatial patterns of lagoon benthic metabolism. Sediment δ 15 N values ranged from 7.2 to 12.1% , and strong spatial gradients were observed in sediment δ 15 N values. The Anselin Local Moran's I identified significant hot spots on the eastern side of BF, and cold spots in upper BSQ ( Figure 1B). δ 15 N values were found to be greater in BF than Brazo de San Quintín (W = 1007, p < 0.01) (Figure 2A), decreased with latitude (τ = −0.45, p < 0.01), and with distance from the ocean (τ = −0.42, p < 0.01) ( Figure 2B).  Figure 3A). Dropping non-significant variables (variable importance for project (VIP) factors <0.8) and choosing a number of factors (8) that minimized RMSE using leaveone-out cross-validation, the predictor variables accounted for 83% of the variability in N stable isotope values. Factors that were highly influential on the regression included Corg/N ratio, residence time, and latitude, whereas sediment S levels were excluded from the regression due to low VIP.  Figure 3A). Dropping non-significant variables (variable importance for project (VIP) factors <0.8) and choosing a number of factors (8) that minimized RMSE using leave-one-out cross-validation, the predictor variables accounted for 83% of the variability in N stable isotope values. Factors that were highly influential on the regression included Corg/N ratio, residence time, and latitude, whereas sediment S levels were excluded from the regression due to low VIP.   Figure 1C). δ 34 S values were found to be greater in BF than Brazo de San Quintín (W = 407, p < 0.01), decreased with latitude (τ = −0.42, p < 0.01), and with distance from the ocean (τ = −0.42, p < 0.01) (Figure 2). Sediment δ 34 S values were modeled successfully as a function of sedimentary and environmental variables ( Figure 3B). Dropping non-significant variables (VIP factors <0.8) and choosing a number of factors (11) that minimized RMSE using leave-one-out cross-validation, the predictor variables accounted for 88% of the variability in S stable isotope values. Factors that were highly influential on the regression included those that scaled with distance from the ocean, including residence time, latitude, and the cost-distance function. The Corg/N ratio was excluded from the regression due to low VIP.
Sediment δ 13 C values ranged from −10.8 to −21.4% , and modest spatial gradients were observed in sediment δ 13 C values. The Anselin Local Moran's I identified significant cold spots in the Arroyo San Simon Delta ( Figure 1D). There was no significant difference in δ 13 C between the two arms of BSQ (W = 585, p = 0.91), although δ 13 C increased with latitude (τ = 0.25, p < 0.01), and with distance from the ocean (τ = 0.26, p < 0.01) (Figure 2). Sediment δ 13 C values were modeled successfully as a function of sedimentary and environmental variables, although the amount of variance explained was substantially lower for δ 13 C than for δ 15 N and δ 34 S ( Figure 3C). Dropping non-significant variables (VIP factors < 0.8) and choosing a number of factors (5) that minimized RMSE using leave-one-out cross-validation, the predictor variables accounted for 55% of the variability in C stable isotope values. Factors that were highly influential on the regression included residence time, % clay, sediment carbonate, and sediment sorting. The Corg/N ratio, depth, d90, and %S were excluded from the regression due to low VIP factors.
Although models that included environmental variables and sediment characteristics were highly predictive of isotopic ratios, it was difficult to identify key drivers-vs. correlates-due to the high level of collinearity among predictor variables (Figure 4). There are essentially two classes of factors in which most variables were strongly positively or strongly negatively correlated. Only two variables-sediment C/N ratio and sediment %P-had modest correlations with most other variables ( Figure 4). As revealed by interpolation maps for the sediment variables measured ( Figures S1 and S2), C/N ratio and sediment P followed different patterns of spatial dependence. For instance, sediment P displayed two minima, at the mid-portions of BSQ and BF, but was low in both the mouth and head locations. The sediment Corg/N ratio followed a similar pattern with high values in mid-BSQ, and lower values near the mouth and head.
Sediment particle size distribution varied from medium silt to fine sand, consistent with a median grain size that ranged from 8 to 190 µm. In general, the coarsest sediment was found near the mouth of the lagoon and in the sub-aqueous portion of the Arroyo San Simon Delta, whereas the finest sediment was found near the heads of BF and Brazo de San Quintín ( Figure S1). Sediments collected near the mouth of the lagoon were moderately well sorted, and fine skewed, grading to very poorly sorted, with symmetrical distributions, near the lagoon heads ( Figure S1). Kurtosis-which refers to the peakedness of sediment distribution curves-was leptokurtic (more peaked) near the mouth of the lagoon, and platykurtic (less peaked) near the lagoon heads ( Figure S1). Plotting sediment mean grain size against sorting revealed that sediments were generally finer with distance from the lagoon mouth, and that sediments further from the mouth were typically more poorly sorted relative to their size than sediments found closer to the mouth of the lagoon (Figure 6). In addition, some sediment samples collected upstream of the tidal restriction indicated deposition under closed basin conditions ( Figure 6).   (Table 1), including seston, Ulva species, Spartina foliosa, Zostera marina, Salicornia pacifica, and Crassostrea gigas.
Sediment particle size distribution varied from medium silt to fine sand, consistent with a median grain size that ranged from 8 to 190 µm. In general, the coarsest sediment was found near the mouth of the lagoon and in the sub-aqueous portion of the Arroyo San Simon Delta, whereas the finest sediment was found near the heads of BF and Brazo de San Quintín ( Figure S1). Sediments collected near the mouth of the lagoon were moderately well sorted, and fine skewed, grading to very poorly sorted, with symmetrical distributions, near the lagoon heads ( Figure S1). Kurtosiswhich refers to the peakedness of sediment distribution curves-was leptokurtic (more peaked) near the mouth of the lagoon, and platykurtic (less peaked) near the lagoon heads ( Figure S1). Plotting sediment mean grain size against sorting revealed that sediments were generally finer with distance  (Table 1), including seston, Ulva species, Spartina foliosa, Zostera marina, Salicornia pacifica, and Crassostrea gigas. the mouth of the lagoon, and platykurtic (less peaked) near the lagoon heads ( Figure S1). Plotting sediment mean grain size against sorting revealed that sediments were generally finer with distance from the lagoon mouth, and that sediments further from the mouth were typically more poorly sorted relative to their size than sediments found closer to the mouth of the lagoon (Figure 6). In addition, some sediment samples collected upstream of the tidal restriction indicated deposition under closed basin conditions ( Figure 6).

Figure 6.
Depositional domain analysis results for lagoon sediments. Sorting and mean values were estimated using the method of moments. Depositional domains based on [71,72].
Supervised classification of seagrass cover indicated that 36% of the lagoon was covered by seagrass in July 2017 (1325 ha), whereas digitization of oyster racks suggested a minimum of 135,000 Supervised classification of seagrass cover indicated that 36% of the lagoon was covered by seagrass in July 2017 (1325 ha), whereas digitization of oyster racks suggested a minimum of 135,000 m of oyster culture racks in the lagoon. Seagrass classification accuracy was 96% based on assessment of 60 polygons. Assessment of sediments collected from plots covered by seagrass (patchy or continuous) suggests that there were no significant differences in δ 15 N or δ 13 C signature (W = 612, p = 0.83; W = 484, p = 0.09) or variance (χ 2 = 0.04, p = 0.84; χ 2 = 0.09, p = 0.76) in relation to seagrass cover ( Figure 7A). There were also no significant differences in δ 34 S signature between plots having or lacking seagrass cover (W = 299, p = 0.48), but there was a trend towards lower δ 34 S values where seagrass was present, and δ 34 S variance was lower in plots with seagrass present (χ 2 = 4.55, p = 0.03). Comparing plots near seagrass meadows (<20 m) vs. far from seagrass meadows (>20 m), we found lower δ 15 N ratios in plots in or adjacent to seagrass beds (W = 427, p = 0.04), but no difference in stable isotope ratios for δ 34 S or δ 13 C (W = 187, p = 0.08; W = 660, p = 0.74). There was no difference in variance between plots distant from and adjacent to seagrass beds for δ 15 N (χ 2 = 1.06, p = 0.15) or δ 13 C (χ 2 = 0.09, p = 0.76), but a lower variance was found for δ 34 S for plots adjacent to seagrass beds (χ 2 = 5.0, p = 0.02).

Figure 7.
Distribution of sediment bulk stable isotope ratios for samples with and without eelgrass cover (A), for samples near (<20 m) and distant from (>20 m) seagrass meadow (B), and for samples near (<100 m) and distant from (>100 m+) oyster culture (C). Statistics represent the results of a Mann-Whitney-Wilcoxon test for magnitude (W) and the Fligner-Kileen test for homogeneity of variance (χ 2 ).

Discussion
BSQ serves as a model for understanding the impact of upwelling on Pacific coastal lagoons in the absence of high anthropogenic impacts. We first interpret our results in light of what is known about spatial patterns of productivity, respiration, and nutrient processing in BSQ. We then address the implications of our results for understanding the impacts of seagrass cover and oyster culture on

Discussion
BSQ serves as a model for understanding the impact of upwelling on Pacific coastal lagoons in the absence of high anthropogenic impacts. We first interpret our results in light of what is known about spatial patterns of productivity, respiration, and nutrient processing in BSQ. We then address the implications of our results for understanding the impacts of seagrass cover and oyster culture on lagoon biogeochemical processing more generally. We conclude by suggesting directions for future research.
BSQ is strongly affected by upwelling of nutrient and pCO 2 -enriched waters, which are also high in 15 N (8-10% ) because of proximity to one of the world's major centers of open-ocean denitrification [47]. Upwelling of nutrient-rich water drives spatial patterns in productivity within the lagoon. As measured by chlorophyll concentrations, productivity is high near the mouth, but declines with distance up BF and Brazo de San Quintín [73]. Dissolved inorganic N (DIN) follows no consistent pattern within the lagoon [74] but may similarly decline with distance from the mouth during upwelling episodes [75]. The lagoon is net heterotrophic, and deeper areas of the estuary are more heterotrophic than seagrass meadows, marshes, and tidal flats, where respiration and production are approximately equal [50]. Based on stoichiometry between O 2 consumption and DIC production in benthic chambers, approximately half of net respiration has been estimated to occur in association with anaerobic metabolism [74]. DIN deficits suggest the lagoon is net denitrifying [49], and spatial patterns in DIN fluxes and macrophyte isotope ratios have suggested the mouth of the estuary sees high denitrification, whereas the arms of the lagoon have a more even balance between N-fixation and denitrification [8,49,76]. Recent analysis of N-fixation has suggested that heterotrophic non-sulfate reducing bacteria had the largest impact on N-fixation at BSQ [77].
Our work found strong spatial patterns in δ 15 N in BSQ, with higher levels in BF, and δ 15 N ratios that declined with distance from the mouth ( Figure 1B). If in this estuary, sediment δ 15 N can be considered a rough indicator of the balance between denitrification, which enriches the 15 N pool through loss of 14 N to the atmosphere, and N-fixation, which introduces light N isotopic signatures (0% ) [78], our results indicate N-fixation is more common in the upper arms of the lagoon, whereas denitrification is more common near the mouth where oceanic imports of nitrate occur. This is consistent with past research into macrophyte isotopic ratios and spatial patterns of N-fixation at BSQ [8,77]. However, our study also identified higher values for sedimentary δ 15 N in BF in comparison with BSQ, which has not been previously reported [8]. The higher sediment 15 N levels in BF may indicate trophic enrichment via oyster biodeposits, or alternatively, enhanced denitrification. Enhanced denitrification may be occurring due to oyster biodeposits or organic matter that act as labile organic matter sources to fuel heterotrophic denitrification [79]. Sediment S isotopic signatures (−33.4 to −12.1% ) did not reflect the S signatures of organic matter sources (−4.3 to +20.1% , excluding seston), suggesting that sedimentary isotopic ratios are indicative of processing. Although previous studies have not measured sulfate reduction in BSQ, the net heterotrophic status of the estuary [49], in combination with high estimates (50%) of anaerobic metabolism [50], suggest that sulfate reduction is common. Fractionation with sulfate reduction can range up to 42% [80], imparting a highly negative signature to reduced S compounds (e.g., pyrite). If sediment δ 34 S can be considered a rough indicator of sulfate reduction [81], spatial patterns in isotopic signatures suggest less sulfate reduction is occurring in BF, and sulfate reduction is more common in the upper arms of the lagoon. Based on spatial patterns in sediment texture ( Figure S1), which grade from sand near the mouth of the lagoon to silts in its upper arms, higher levels of sulfate reduction might be expected in the upper bays due to lower oxygen permeability. Additionally, higher carbon densities of sediments in the upper bays ( Figure S2) may indicate higher availability of organic matter for sulfate reducers, leading to increased sulfate reduction as reflected in more negative δ 34 S. The "hot spot" of higher δ 34 S values in areas of intense aquaculture may be due to high nitrate abundances that favor denitrification over sulfate reduction as a mode of anaerobic metabolism. Furthermore, the low δ 15 N values coupled with low δ 34 S values in upper BSQ suggest that N-fixation and sulfate reduction may be occurring as a coupled process [18], although recent work in the lower part of the lagoon found that N-fixation was linked to non-sulfate reducing bacteria [77].
Stable isotopes are often used to trace seagrass organic matter in food webs, due to the distinct combination of unusually heavy C isotopic signatures (indicative of bicarbonate assimilation) [82] combined with light S isotope signatures indicative of their rooting in reduced environments [83]. Despite the dominance of seagrass in this estuary (36% of areal extent), stable C isotopes in sediments in BSQ (mean = −18.3% ) were much more similar to presumed phytoplankton values (−20% , [84] and measured seston samples, mean = −19.7% ) rather than those of Zostera marina (mean = −11.2% ). In addition to organic matter source, stable C isotopes are also known to reflect interlinking processes of productivity, depth, or C source (bicarbonate vs. CO 2 ) in coastal marine settings [3]. Under high productivity, partial CO 2 limitation results in lower discrimination against 13 C and higher δ 13 C signatures [85]. In the coastal ocean, light is often production-limiting for phytoplankton, such that shallower areas have less negative δ 13 C values more typically associated with high production, and deeper areas have lighter δ 13 C signatures. Finally, there are strong gradients in pCO 2 in BSQ, and isotopic signatures of δ 13 C and HCO 3 are distinct [86]. Because we see lighter C isotopic signatures in areas of high productivity ( Figure 1D) (rather than the opposite) and see no significant relationship between δ 13 C and depth (Figure 4), we can rule out these factors as providing important controls on spatial patterns in δ 13 C. We also see no consistent relationship between macrophyte δ 13 C values and distance to the lagoon mouth (τ = −0.06, p = 0.57), which we might expect if spatial gradients in pCO 2 or HCO 3 played a large role in determining sediment δ 13 C values.
In addition to the specifics of macronutrient cycling in BSQ, this work has broader implications for understanding the role of seagrass presence and oyster culture in coastal lagoons. Seagrasses have often been associated with enhanced levels of sediment N fixation [87], because both their exudates and the organic matter trapped by seagrass meadows serve as fuel for heterotrophic diazotrophs [84,88].
Although seagrass has been associated with N-fixation in BSQ previously [77], our study did not identify associations between seagrass cover with the depressed values of δ 15 N that would result from enhanced N fixation at either the plot or landscape scale (Figure 7). Seagrasses may also enhance sulfate reduction through exudates [89] or decrease sulfate reduction through root zone oxidation [90]. Although we found no evidence for lower δ 34 S values in sediments collected in or nearby seagrass, the range of sediment S isotopic values was consistently reduced in seagrass covered compared to in bare sediment ( Figure 7A). Finally, while the seagrasses we measured tended to have C and S isotopic signatures that were distinctive based on high δ 13 C and low δ 34 S, the S isotopic signatures were altered by processing and not preserved in sediments ( Figure 1C), rending these values not useful for tracing seagrass organic matter in seagrass meadows.
Another factor thought to shape coastal nutrient cycling is the establishment of shellfish aquaculture. In this study, we found consistently higher δ 15 N values near areas of oyster culture. Previous research has found strong relationships between oyster culture and denitrification [91], and shell material alone has been observed to enhance denitrification [92]. Gene-based analysis of the microbiome of Crassostrea virginica found denitrifying bacteria in the oyster tissue and shell, further supporting oyster reefs as hotspots for denitrification [93]. No significant difference was found in S isotopes from oyster culture areas in comparison with non-culture areas ( Figure 1C). This pattern suggests that oyster culture may not be enhancing sulfate reduction, which has been noted previously in association with Crassostrea gigas culture [94], and can be a matter of concern in terms of estuarine water quality [95].

Conclusions
This study has demonstrated the role that spatial analysis of sediment stable isotope data can play in identifying potential patterns in lagoon metabolism related to N and S processing, and the dispersion of organic matter. Coastal ocean isoscape mapping has proved useful in identifying spatial patterns in productivity, DIC dynamics, and nutrient pollution [6][7][8]. However, as a tool it is not well exploited in comparison with more standard approaches, such as water quality analysis, or even nutrient processing measures, which are substantially more limited in terms of the time window they reflect. This research suggests that stable isotope mapping can provide meaningful insight into the assessment of estuarine metabolism, or can act as a stand-alone tool for generating hypotheses, identifying the influence of spatial gradients, or suggesting prime locations for investigation of microbial abundance or function.