Seasonal patterns of benthic−pelagic coupling in oyster habitats

Oysters enhance benthic−pelagic coupling in coastal systems by moving large quantities of suspended particulates to the sediments, stimulating biogeochemical processes. Recent research efforts have focused on quantifying the impact of oysters on coastal biogeochemical cycling, yet there is little consensus on how oysters influence processes across systems. A potential driver of this variance is availability of organic material suspended in the water column and subsequent loading to sediment by oysters. Here, we measured fluxes of sediment di-nitrogen (N2-N), ammonium (NH4), combined nitrate-nitrite (NOx), and phosphate (PO4) in spring, summer, and fall at 2 oyster reefs and 1 farm in a temperate estuary (Narragansett Bay, Rhode Island). We then linked these fluxes with patterns of water column primary production. Nitrogen removal from the system was highest in spring, when we detected net sediment denitrification (48.8 μmol N2-N m−2 h−1) following a winter−spring diatom bloom. In contrast, we measured sediment N2 fixation in fall (−44.8 μmol N2-N m−2 h−1) at rates nearly equivalent to spring denitrification. In the summer, we measured a nearly net zero sediment N2-N flux (−2.7 μmol N2-N m−2 h−1). Recycling of nitrogen to the water column was consistent across seasons, composed almost exclusively of NH4. These results demonstrate that sediment nitrogen cycling in oyster habitats is dynamic and can change rapidly based on seasonal patterns of productivity. At carrying capacity, the impact of oysters on nitrogen cycling is large and should be considered during efforts to increase oyster populations through aquaculture or reef restoration.


INTRODUCTION
Large oyster populations characterized many estuaries on the Atlantic Coast of North America prior to European colonization (Nixon 1997, Beck et al. 2011, Zu Ermgassen et al. 2013, Rick et al. 2016. Globally, oyster reefs have declined by approximately 85% of their historic extent (Beck et al. 2011), and less than 6% of oyster reefs on the Northeastern Coast of the United States remain (Zu Ermgassen et al. 2012). Increasingly, there are efforts to restore oyster reefs and develop the oyster aquaculture industry to recover both economic and ecological benefits. Such ecological services include provision of habitat for other marine species, protection from storm surge, and prevention and reduction of the symptoms of eutrophication (Jackson 2001, Newell et al. 2005, Bernhardt & Leslie 2013. Additionally, oysters play an important role in regulating benthic−pelagic coupling processes in coastal ecosystems (Newell et al. 2005, Smyth et al. 2013a. Specifically, oysters enhance rates of organic matter (OM) deposition to the benthos, thereby increasing dissolved nitrogen (N) and phosphorus (P) recycling in coastal systems (Boucher & Boucher-Rodoni 1988, Dame et al. 1992. Oysters not only enhance elemental cycling rates in coastal systems, but can also change the pathways dominating these rates. For example, oysters may stimulate denitrification, a process that converts bioavailable nitrate (NO 3 − ) and nitrite (NO 2 − ; collectively NO x ) into inert di-nitrogen (N 2 ) gas.
In coastal systems polluted with anthropogenic N, denitrification can serve to remove up to 50% of excess N from the system (Seitzinger 1988), and there is substantial interest in the potential for oyster restoration projects and aquaculture to promote denitrification (Newell et al. 2005). Recent efforts to quantify how rates of sediment denitrification change in the presence of oysters yield conflicting results. Some studies report enhanced denitrification in sediments beneath oysters in reef and aquaculture systems (Kellogg et al. 2013, Smyth et al. 2015, Lunstrum et al. 2018, while others have found no difference (Higgins et al. 2013, Testa et al. 2015, Westbrook et al. 2019. Similarly, there is large variance in the amount by which denitrification is enhanced when higher rates are found in oyster habitats, with some studies reporting enhancement of less than 100 μmol N 2 -N m −2 h −1 (Hoellein et al. 2015, Lunstrum et al. 2018, others several hundred μmol N 2 -N m −2 h −1 , Smyth et al. 2016, Ray et al. 2020, and even reports of en hancement over 1000 μmol N 2 -N m −2 h −1 (Kellogg et al. 2013). Perhaps the reason for this disparity is the assumption that oysters will enhance sediment denitrification similarly regardless of the properties of the system they are in (e.g. dissolved nutrient concentrations, temperature, turbidity, etc.). Another reason for this disparity could rest in the complexity of the N cycle.
Coastal sediment N cycling is dynamic, with multiple, often competing, processes co-occurring (Fulweiler et al. 2013, Kuypers et al. 2018, Foster & Fulweiler 2019, Ray et al. 2020. For example, coastal sediments can rapidly (days to months) switch from N 2 fixation -which introduces new N to the system -to denitrification following the addition of labile OM (Fulweiler et al. 2007(Fulweiler et al. , 2008(Fulweiler et al. , 2013. As this OM is used up and only recalcitrant OM remains, net N 2 fluxes decline (Fulweiler et al. 2013) and may become negative as N 2 fixation proceeds more rapidly than denitrification (Fulweiler et al. 2008). In highly productive systems and/or seasons, we anticipate that denitrification will be promoted in oyster habitats as oysters move labile OM to the sediments below in biodeposits consisting of both feces and pseudofeces. In less productive systems and/or seasons, oysters may be more likely to consume and digest any suspended particles that are processed during filter-feeding, yielding fewer biodeposits, thus reducing OM loading to the sediments beneath, favoring N 2 fixation. We hypothesize that net sedi-ment N 2 fluxes in oyster habitats track primary production in the water column, with net denitrification following large blooms, followed by a rapid switch to N 2 fixation when productivity and availability of labile OM added to the sediments via oyster biodeposits declines (Burgin & Hamilton 2007, Fulweiler et al. 2013, Hardison et al. 2015. Sediment O 2 consumption and phosphate (PO 4 3− ) release can be enhanced in oyster habitats (Mazouni et al. 1996, Kellogg et al. 2013 and are also likely to follow patterns of water column productivity. Greater OM loading to sediments can drive aerobic decomposition, reducing O 2 concentrations and leading to PO 4 3− release (Ingall & Jahnke 1994, Foster & Fulweiler 2019). Thus, we expect greater O 2 consumption and PO 4 3− release immediately following blooms.
Quantifying seasonal patterns of benthic−pelagic coupling in oyster habitats will help us understand biogeochemical cycling in coastal systems prior to over-exploitation of oyster populations and allow for smart decision making as oyster populations are enhanced through restoration and development of aquaculture. Here, we took advantage of seasonal patterns of primary production in Narragansett Bay to investigate how these patterns regulate sediment biogeochemical processes in oyster habitats. We assessed whether seasonal patterns of sediment biogeochemical processes followed bay-wide patterns of primary production. We further tested whether the magnitude of these biogeochemical processes could be predicted by water column and sediment chlorophyll a (chl a) and pheophytin (pheo) concentrations, where chl a indicates the relative abundance of live phytoplankton, and pheo indicates the relative abundance of degraded chl a (Yentsch & Menzel 1963), indicative of dead or dying phytoplankton. To accomplish this goal, we collected sediment and water column chl a and pheo samples to compare with a long-term weekly data set of sur face pigment concentrations (The Narragansett Bay Long-Term Plankton Time Series, https:// web. uri. edu/ gso/ research/ plankton/), and incubated sediment cores to quantify rates of sediment O 2 consumption, and N and PO 4 3− fluxes across the sediment−water interface.

Site description and sampling scheme
Narragansett Bay is a shallow (mean depth: 8.6 m) temperate estuary. The Bay runs from north to south and drains to Rhode Island Sound. The Providence River in the north is the main freshwater source, but daily freshwater inflow is low (annual mean: 100 m 3 s −1 ; Pilson 1985). The average water retention time is approximately 26 d (Pilson 1985). Concentrations of dissolved N and P are typically highest near points of freshwater inflow and decrease along a north to south gradient (Smayda & Borkman 2008). Weekly measurements of water column chl a and pheo have been collected at 41.567°N, 71.384°W since the 1950s, and data collected since 1999 is available online as part of The Narragansett Bay Long-Term Plankton Time Series (hereafter Plankton Time Series). Phytoplankton blooms in the Bay tend to be synchronous bay-wide events, regardless of the point of initiation (Smayda & Borkman 2008). Primary production in the system has historically been characterized by large winter−spring blooms, but the period of peak bloom has become more variable over the last 50 yr ).
Like other coastal systems near large human populations, Narragansett Bay has experienced significant change over the past 400 yr. Prior to European colonization, oysters covered a large portion of the Bay. Following over-harvest and mortality associated with disease, pollution, and siltation, less than 1% of natural populations remained at the start of the 21 st century (Zu Ermgassen et al. 2012, Schumann 2015. Modeling estimates suggest the Bay could currently sustain 625 times more oysters before reaching ecological carrying capacity (Byron et al. 2011). More recently, there have been efforts to increase oyster populations within Narragansett Bay through restoration and the expansion of the oyster aquaculture industry (Beutel 2017). Between 2011 and 2016, sales of cultured oysters raised in the state of Rhode Island (which includes Narragansett Bay and several smaller coastal lagoons that are at or near the legal limit for aquaculture development) have increased by 92% (Beutel 2017). Beginning in 2006, wastewater treatment plants in the watersheds that drain to Narragansett Bay implemented ad vanced N removal technologies. By 2012, the N load was 50% less than in 2006, and in 2016 the total N load was 66% less than in 2006 (Narragansett Bay Estuary Program 2017, Oczkowski et al. 2018). During this same time frame, Narragansett Bay has also shown symptoms of a system responding to a changing global climate (Fulweiler et al. 2007), with a reduced magnitude of the winter−spring diatom bloom, which has only occasionally occurred in recent years. Warmer winters in the region are the main driver of a reduced winter− spring bloom, as warm winters appear to enhance rates of grazing on phytoplankton, increase the number of cloudy days which subsequently reduces the amount of light available for photosynthesis, and make the water column more stable, leading to less mixing and nutrient availability in the upper portion of the water column , Fulweiler et al. 2013, Fulweiler & Heiss 2014.
We collected sediment and water column samples from 3 oyster habitats in Narragansett Bay between July 2016 and October 2017 ( Fig. 1). Two sampling sites were on the western side of the bay. Allen Harbor is an oyster farm located directly on the main channel of the Bay that employs the rack and bag method of culture, where oysters are suspended in mesh bags approximately 30 cm from the sediment surface. Bissel Cove is a small embayment that contains both natural and restored oyster reefs, with restoration taking place within the 5 yr prior to our sampling. The third site, Town Pond, is located in the northeastern portion of Narragansett Bay in Mt. Hope Bay. The embayment contains several small restored oyster reefs which were constructed between 2008 and 2014. The embayment is surrounded by saltmarsh and is in close proximity to a golf course. At both reef sites, oysters were in direct contact with the sediments.

Water column chl a sample collection and analysis
We collected water column chl a and pheo samples monthly at each site during spring (April, May), summer (June, July, August), and fall (September, October). Each month between July 2016 and October 2017, we collected samples every 20 min over a 3 h period, and once a season, we collected samples every hour for 8 h. During each sampling event, we collected duplicate samples by filtering site water through GF/F (0.7 μm) filters, which were stored in a cooler on ice until return to the lab where they were stored at −80°C until analysis. Samples were extracted in 10 ml of 90% acetone and kept at 4°C for 24 h until fluorometric analysis (Arar & Collins 1997). Immediately following fluorometric analysis, we added 25 μl of 0.30 HCl to each sample and then waited 3 min before repeating the analysis in order to estimate pheo concentration (Arar & Collins 1997).

Sediment collection and flux measurement
We hand-collected triplicate sediment cores from each site in summer (August) and fall (October) 2016, and spring (April) and summer (July) 2017 using clear PVC tubes that were 28 cm in length and had a cross sectional area of 0.00785 m 2 . Approximately 12 cm of sediment was collected. Cores were collected from patches within the reef where no oysters were present (Bissel Cove and Town Pond) and from directly beneath oyster cages (Allen Harbor). Sediments in Bissel Cove and Town pond were silty and mucky and made up of fine particles, while sediments at Allen Harbor were made up of fine sand and silt. Immediately after collection, we added site water to the top of each core before capping them, placing them in a cooler, and transporting them to the laboratory for incubations and analysis. We also collected unfiltered site water in 20 l acid-washed carboys for use during the incubations. Upon return to the lab, we moved the cores into a water bath in an environmental chamber, both maintained at the in situ temperature of the field site. We removed the caps and gently bubbled the cores overnight. Prior to beginning incubations the next day, we siphoned off the water overlaying the sediments in each core and replaced it with water from the carboys. We also filled an empty core with unfiltered site water to be incubated alongside the sediment cores to account for any fluxes occurring in the water column. In cases where there was a flux in the water control core, we subtracted this flux equally from the fluxes measured in the sediment cores. This adjustment was made for 14 of 42 cores for N 2 -N fluxes, 30 of 42 cores for O 2 flux, and for all cores for ammonium (NH 4 + ), NO x , and PO 4 3− fluxes. Immediately prior to capping each core and beginning incubations, we measured salinity using a conductivity probe (CDC401; Hach) and dissolved oxygen (DO) using an optical sensor (LDO101; Hach). We then removed 30 ml of water using an acidwashed plastic syringe and filtered this sample through a Whatman GF/F filter (0.7 μm) into an acidwashed polyethylene bottle, which was immediately frozen until later analysis for dissolved inorganic N (DIN) and PO 4 3− . Each core was then capped with a clear acrylic lid equipped with a magnetic stir bar to keep water overlying the sediments gently mixed (40 rpm), an inflow port connected to one of the sitewater carboys, and an outflow port for sample collection (Ray et al. 2019a).
We collected duplicate water samples for analysis of N 2 -N and O 2 concentration at 5 time points, which were spaced to allow DO in each core (except the site-water control) to drop at least 2.0 mg l −1 over the course of the incubation without becoming hypoxic (concentration below 2.0 mg l −1 ) (Foster & Fulweiler 2014, Ray et al. 2019b). Incubations typically lasted 5−10 h, with each sample collected 1−2 h apart. At each sample collection point, we filled 12 ml glass Labco exetainers (n = 2) from the bottom, allowing the exetainers to overflow 3 times the full volume. We then added 25 μl of saturated zinc chloride solution to each in order to stop any biological activity and preserve the sample, then capped them with gas tight cap septa. As samples were collected, both the inflow and outflow ports to the core were opened so that site-water from the carboy could replace the volume of water collected in the sample. Following the last time point, we removed the core lid and again collected samples for DIN and PO 4 3− analysis. Total water replacement in each core over the course of the incubation was less than 15%. Lights remained off for the duration of the incubations, except for short periods to collect samples, check DO levels in the cores, and to check for bubble formation. We recorded a positive DO flux from one core collected from Town Pond in summer 2017, indicating there was a leak allowing for gas exchange with the incubation chamber. We excluded all fluxes from this core from our analyses.
The day after the incubation, we siphoned the water from each core in order to collect samples for analysis of sediment properties. Using an acid-washed 60 ml plastic syringe with the tip cut off, we collected the top cm of sediment and stored these samples in acid-washed centrifuge tubes at −20°C until analysis for porosity, density, % OM, and C:N. We repeated this process using a deionized water (DI)-rinsed syringe and DI-rinsed centrifuge tubes to collect samples for analysis of sediment chl a and pheo content. These samples were stored at −80°C until analysis.

Sample analysis and sediment flux calculation
We quantified sediment porosity and bulk density using standard methods (Nielsen et al. 2000). Sediment OM was measured as the percent of mass lost on ignition at 500°C, and sediment total C and N and C:N were quantified using a Micro Vario Elemental Analyzer (Elementar Americas). We quantified sediment chl a concentration by extracting a known mass of sediment in 92% acetone, then sonicated the sample for 30 s and stored it at 4°C for 24 h until fluorometric analysis (Arar & Collins 1997, Fagherazzi et al. 2014. Sediment pheo was estimated by adding 25 μl of 0.300 HCl to the sample immediately following fluorometric analysis, waiting 3 min, and then repeating the fluorometric analysis.
We quantified net N 2 -N and O 2 fluxes using a quadrupole membrane inlet mass spectrometer (MIMS; Bay Instruments) following the N 2 :Ar method (Kana et al. 1994). N 2 -N concentrations in each sample were calculated by multiplying the sample N 2 :Ar ratio determined by the MIMS by the theoretical Ar concentration in the sample given its temperature and salinity (Weiss 1970, Colt 1984. We estimated sample O 2 concentrations in the same way, but instead used O 2 :Ar determined by the MIMS.
DIN and dissolved inorganic P (DIP) samples were analyzed using high-resolution digital colorimetry (Seal Autoanalyzer 3; SEAL Analytical). Minimum detection limits during sample analysis were 0.80 μM for NH 4 + , 0.013 μM for NO x , and 0.010 μM for PO 4 3− . We estimated fluxes of N 2 -N and O 2 using a linear regression of gas concentration over time (h). When the regression calculation R 2 ≥ 0.65 and p ≤ 0.10, we multiplied the slope of the regression line by the core volume, and then divided by the core cross sectional area to estimate the flux in μmol m −2 h −1 . If R 2 ≤ 0.65 or p ≥ 0.10, the flux was assigned a value of zero (Prairie 1996, Foster & Fulweiler 2016. For DIN and DIP, we calculated the flux as the difference between the concentrations in the initial and final samples, scaled by the core volume.

Statistical analysis and data availability
We conducted all statistical analyses and made all figures using R statistical software version 3.6.1 (R Core Team 2014) and the 'lme4' (Bates et al. 2015), 'fitdistrplus' (Delignette-Muller & Dutang 2015), 'emmeans' (Lenth 2018), 'ggplot2' (Wickham 2016), 'dplyr' (Wickham et al. 2019), and 'plyr' packages (Wickham 2011). The map of site locations was made using the 'ggspatial' (Dunnington 2018) and 'sf' (Pebesma 2018) packages and the NOAA Medium Resolution Shoreline shapefile (https:// shoreline. noaa. gov/ data/ datasheets/ medres. html). We considered the results of all statistical tests to be significant when p ≤ 0.05, and report all rates and concentrations as the mean ± SE. In addition to performing analyses on fluxes of individual compounds, we also analyzed the sum of all measured N fluxes (ΣN), which we calculated by adding the N 2 -N, NH 4 + -N, and NO x -N fluxes measured for an individual core. Finally, we considered each sediment core and the fluxes and properties associated with that core to be independent samples when assessing predictive variables and relationships between sediment properties and fluxes.
We tested whether chl a and pheo concentrations in our hand-collected samples could be predicted by values reported in the Plankton Time Series using a linear regression approach. First, we calculated the mean daily value for chl a and pheo for each sampling date. Next, we matched those values with values from the nearest date in the Plankton Time Series. If our hand sampling fell directly between 2 time series dates, we used the average of the 2 values in the time series. To test for correlation between sediment chemical properties, we used Pearson correlations.
To compare fluxes between seasons, we used a mixed model approach (Zuur et al. 2009). To begin, we determined whether our measured ΣN, N 2 , NH 4 + , NO x , PO 4 3− , and O 2 fluxes best fit normal, log normal, or gamma distributions using the 'fitdistrplus' package. We found that N 2 -N, NH 4 + , and PO 4 3− fluxes were best described by normal distributions, ΣN and O 2 fluxes best fit gamma distributions, and NO x fluxes were log normal. In order to best meet the assumptions of mixed models, we shifted the ΣN data by adding one plus the lowest flux value to every flux, and multiplied the O 2 fluxes by −1, so that all flux values were positive. For the NO x fluxes, we shifted the data so all values were positive in the same way as for the ΣN data before performing a square root transformation. We then compared flux -es between seasons by creating generalized linear mixed models (GLMMs), using season as a fixed effect and sampling site as a random effect, then used least square means tests to conduct pairwise comparisons.
To investigate whether water column and sediment chl a and pheo concentrations could predict sediment fluxes, we used a linear regression approach. For the water column regressions, we used the reported water column chl a, pheo, or sum of chl a and pheo concentrations from the Plankton Time Series the week before the sediment incubation, the mean from the previous month, and the mean from the previous 3 mo as the independent variable, and the mean flux for a given analyte from that incubation as the dependent variable. We considered the best model for each flux to be the one with the highest R 2 value. For the sediment regressions, we used the chl a and pheo concentrations as quantified from samples collected immediately following incubations, and those concentrations specific to each individual sediment core.

Seasonal water column and sediment properties
Dissolved NH 4 + and NO x concentrations were relatively low year round, with the highest concentrations in the spring, specifically at Town Pond (Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m652 p095 _ supp. pdf). Water column chl a concentrations followed the classic temperate pattern of a winter−spring bloom with chl a peaks in January and February and a secondary bloom in July and August in 2016 and 2017 (Fig. 2).
While the concentrations of hand-collected chl a samples at each site tended to follow those recorded as part of the long-term time series (Fig. 2), we only found one statistically significant linear relationship between the pigment concentration of handcollected samples relative to concentrations reported in the long-term time series (Table 1). Sediment density and porosity were similar across seasons (Table S2). Across all sites, the amount of OM was highest in spring (3.97 ± 1.00% OM), followed by fall (3.65 ± 1.22% OM), and lowest during summer sampling (2.74 ± 0.57% OM; Table S2). The   Table S2). Sediment chl a concentrations followed the same pattern as sediment C:N, decreasing from spring to fall (spring: 3857.9 ± 642.6 μg chl a m −2 , summer: 2828.1 ± 585.4 μg chl a m −2 , fall: 2367.9 ± 562.4 μg chl a m −2 ; Table S2). Sediment pheo patterns were similar, with highest concentrations in spring (5650.7 ± 1202.9 μg pheo m −2 ), a marked drop in summer (2781.1 ± 559.1 μg pheo m −2 ), and then a slight increase in fall (3770.4 ± 1007.5 μg pheo m −2 ; Table S2). Many sediment properties were significantly correlated with one another (Table 2). Pheo content in sediments was significantly related to the total amount of OM, % C, % N, and C:N. We found no significant relationships between chl a and any other sediment property.

Chl a and pheo concentrations as predictors of biogeochemical processes in oyster habitats
There were few statistically significant relationships between mean water column chl a and pheo concentrations and sediment fluxes (Table S3). There were no significant relationships between mean chl a and pheo concentration and sediment O 2 , PO 4 3− , N 2 -N, or NO x fluxes. However, ΣN fluxes could be predicted by the mean water column pheo concentration from the previous month (p = 0.01, R 2 = 0.51), with lower ΣN fluxes at high pheo concentrations. Sediment NH 4 + fluxes were best predicted by the mean of the sum of chl a and pheo from the previous month (p = 0.04, R 2 = 0.30), again with high concentrations associated with lower NH 4 + fluxes. Sediment chl a and % OM were not significant predictors of any fluxes (Table S4). However, sediment pheo concentration was significantly related to sediment ΣN (p = 0.03, R 2 = 0.12) and N 2 -N (p < 0.01, R 2 = 0.21) fluxes, but not NH 4 + or NO x (Fig. 4, Table S4).

Seasonal patterns of benthic−pelagic coupling in oyster habitats
Sediment biogeochemical processes in oyster habitats followed seasonal patterns of water column productivity, with highest ΣN fluxes in April (spring; Fig. 3 August (summer) as sediment metabolism increased with temperature, before dropping dramatically in October (fall; Fig. 3) as temperature fell and less labile OM was available. Sediment PO 4 3− release was highest in the summer when O 2 demand was greatest.
Mean sediment NH 4 + release was not significantly different across seasons (Fig. 3). In Narragansett Bay, most N 2 production is from coupled nitrification− denitrification (Seitzinger et al. 1984, Seitzinger & Nixon 1985, Fulweiler & Nixon 2012. It is likely that a significant portion of the NH 4 + produced via remineralization in the spring was quickly used in coupled nitrification−denitrification, reducing total net NH 4 + efflux. As the availability of labile OM that could be used during denitrification decreased during summer, NH 4 + production fell and coupled nitrification− denitrification simultaneously declined. Net N 2 -N fluxes were similar between sites in spring and fall, but more variable in summer, with some differences between sites (Fig. 3). For example, in summer 2016, we measured only net zero fluxes at Allen Harbor and Bissel Cove, and denitrification in Town Pond, with fluxes measured during the start of a bloom (Figs. 2 & 3). In summer 2017, fluxes were measured before the secondary summer bloom, and both Town Pond and Allen Harbor displayed net N 2 fixation, while Bissel Cove displayed denitrification (Fig. 2 & 3).
We also considered sediment N fluxes in the context of ΣN. Using ΣN allowed us to consider the influence of oyster habitats on the overall sediment N metabolism and amount of N that is processed, regardless of its form. Seasonal patterns of ΣN fluxes were driven by N 2 -N fluxes, with N 2 -N efflux (deni- , and NO x fluxes at oyster habitats in Narragansett Bay. Distinct letters indicate seasons significantly different (p ≤ 0.05) from one another following least square means tests. Each point indicates a single flux measurement, with those above the zero-line indicating release to the water column, and those below indicating sediment uptake. The solid line in the middle of each box is the median, the lower and upper edges of the boxes indicate the 25th and 75th percentile, and the end of each whisker indicates the smallest and largest value within 1.5 times the interquartile range trification) in the spring, a net zero N 2 -N flux in the summer, and N 2 -N consumption (N 2 fixation) in the fall (Fig. 3). As we hypothesized, denitrification was highest in the spring following biodeposition after (or from) the winter−spring bloom. We found that sediment pheo concentration was a significant predictor of net N 2 -N fluxes, with higher N 2 -N release at greater pheo concentrations (Fig. 4). Denitrification tended to be high when labile OM was available, as in the spring following high biodeposition during the winter−spring diatom bloom. As the labile OM in sediments is consumed, leaving behind recalcitrant OM, denitrification slows and N 2 fixation should dominate. Supporting this idea, Fulweiler et al. (2008) observed a 'N 2 fixation threshold' in Narragansett Bay, where sediments switched to N 2 fixation when OM deposition fell below 0.3 g C m −2 d −1 . We did not measure rates of C deposition, but calculated an estimate. Riemann et al. (1989) estimated a g C:g chl a ratio of between 27 and 67 for mixed phytoplankton communities with various growth rates. Using this ratio, we can estimate summertime water column C the month before sampling as 244.6− 607.2 μg C l −1 in summer 2016 and 184.7−458.3 μg C l −1 in summer 2017. Using an oyster filtration maximum of 0.17 m 3 g dry weight [DW] -1 d −1 (Ehrich & Harris 2015), and an oyster density of 5.93 g DW m −2 (Byron et al. 2011), we can estimate that approximately 1.01 m 3 d −1 of water are filtered per every m 2 of oyster habitat. If the oysters consumed no C in that water, moving all of it instead to the sediments (an unlikely scenario), C loading to sediments was approximately 0.25−0.61 g C m −2 d −1 in 2016 and 0.19−0.46 g C m −2 d −1 in 2017. In both cases, the lower end estimate of the C:chl a conversion is below the N 2 fixation threshold of 0.3 g C m −2 d −1 before even including oyster metabolic C re quirements. Newell (1988) estimated that oysters have a summertime metabolic requirement of 17.3 mg C g DW -1 d −1 . If there are 5.93 g DW m −2 , this equates to 0.1 g C m −2 d −1 , further reducing the summertime sediment C load. We recorded both net denitrification and net N 2 fixation during our summer measurements, reflective of the high and low ends of our sediment C loading estimates, and capturing natural variation around this important threshold.
We expected to record denitrification in the fall following the secondary summer bloom, but were surprised to measure such strong N 2 fixation, especially considering our flux measurements were made within a few weeks of the bloom (Figs. 2 & 3). There are a few ways this pattern may be explained in a Fig. 4. Relationship between sediment pheophytin concentration (top 0−1 cm) and fluxes of nitrogen compounds between the sediment and water column. R 2 values indicate the result of linear regressions. Each point represents an individual measurement, with values above zero indicating release from the sediment to the water column, and points below indicating sediment uptake. Shaded region around the best fit line: 95% CI biogeochemical context, taking seasonal patterns of sediment physical−chemical properties and oyster metabolism into account. Despite a similar % OM content in spring and fall measurements, sediment C:N ratio was lowest in spring and highest in fall, indicating that OM added to sediments in spring was more labile than in fall. This could be driven by differences in the phytoplankton community and potential OM quality, or high rates of microbial activity in summer decomposing OM more efficiently and leaving behind only recalcitrant OM.

Water column pigments do not predict sediment fluxes in oyster habitats
Decades of data collection in Narragansett Bay show a tight coupling between mean summer chl a and sediment N 2 -N flux (e.g. Fulweiler & Heiss 2014), and we were surprised not to record similar patterns here (Table S4). Fulweiler & Heiss (2014) found a strong predictive relationship (p < 0.01, R 2 = 0.87) between summer water column chl a concentration and sediment N 2 -N fluxes in the main stem of Narragansett Bay. Using their predictive equation (net N 2 -N flux = 18.3 × chl a concentration − 33.4), we should have measured a net N 2 -N flux of 114.5 μmol m −2 h −1 in summer 2016 but instead recorded a rate about a quarter of this (35.05 μmol m −2 h −1 ). In summer 2017, the predicted flux was 125.81 μmol m −2 h −1 , but we measured net N 2 fixation (−60.28 μmol m −2 h −1 ). There are several factors that may explain the difference in the patterns found here compared to the literature. First, our hand-collected water column pigment concentrations did not closely match with those of the long-term Plankton Time Series used by Fulweiler & Heiss (2014). Second, oysters may have consumed significant portions of the chl a and pheo from water column particulates before moving biodeposits to the sediment in the summer and fall, changing the loading or quality of deposited OM as described previously. Third, it is possible that the dominant phytoplankton species varied across seasons, and oysters selectively ingested or rejected particulates of different quality OM. Finally, our measurements were made in shallow water (1−2 m), and those from Fulweiler & Heiss (2014) were made at depth (~7 m) in the main stem of the Bay.
Our measurements of water column chl a and pheo from embayments did not match those from the Narragansett Bay Time Series, which is collected in the main stem of the Bay. At all of our sampling sites, it is likely that consumption of phytoplankton by oysters played a role in this difference. On almost every hand-sampling event, the concentration of handcollected samples was lower than time series values (Fig. 2). Two of our sampling sites (Bissel Cove and Town Pond) are connected to the main stem of the Bay by small channels, so any phytoplankton that entered the embayment on a rising tide was likely quickly consumed and was depleted by the time the tide fell. At the third site (Allen Harbor), a series of oyster cages separated the main stem of the Bay from our chl a collection site, and it is possible these oysters consumed phytoplankton from the water column before it reached our sampling site, reducing chl a and pheo concentrations. Regardless, using the time series provides an overview of seasonal patterns of phytoplankton abundance and can still be used to help us understand seasonal dynamics of benthic− pelagic coupling in oyster habitats within the system.
Our measurements of % OM addresses the total amount of OM available in sediments. It is likely that as phytoplankton and detritus pass through the oyster digestive tract, the most labile OM is digested first, leaving behind recalcitrant OM. In fact, there is relatively high revival capacity for diatoms that have passed through an oyster (Barille & Cognie 2000). As described in the previous section, sediment C loading was near the N 2 fixation threshold, and denitrification requires labile OM, while N 2 fixation can proceed with more recalcitrant OM. If oysters remove the labile portion of food particles and deposit the recalcitrant portion, N 2 fixation may likely still dominate, even if C loading is above the N 2 fixation threshold. Finally, variations in the phytoplankton community, and the associated OM quality, oyster digestive preferences, and sediment loading may also vary seasonally.

Comparison with studies in other coastal systems
The N 2 -N flux rates we measured in summer (−2.7 μmol N 2 -N m −2 h −1 ) are much lower than summertime rates measured beneath aquaculture (553.6 μmol N 2 -N m −2 h −1 ; Ray et al. 2020) in nearby Ninigret Pond (< 50 km distance), and lower than nearly all other studies that report N 2 -N fluxes in oyster habitats. Our measured rates of denitrification in spring (48.8 μmol N 2 -N m −2 h −1 ), however, are within the range reported in other studies (Higgins et al. 2013, Hoellein et al. 2015, Onorevole et al. 2018. The high rate of N 2 fixation we found in fall (−44.8 μmol N 2 -N m −2 h −1 ) is the highest mean N 2 fixation rate ever reported for oyster habitats. We are not the first to report N 2 fixation in oyster habitats (Mortazavi et al. 2015, but are the first to find it consistently. The majority of studies measuring net sediment N 2 fluxes in oyster habitats have been conducted in a few coastal systems on the Atlantic Coast of the USA. Many take place in the summer alone, and often for 1 or 2 yr. This limited temporal sampling may drive the large variance in reported rates of denitrification in these studies, due to seasonal and annual patterns of OM availability in the water column. We can compare the work conducted here with seasonally reported rates from the Chesapeake Bay and Bogue Sound -2 systems with multiple studies that report seasonal rates of sediment denitrification in oyster habitats. Kellogg et al. (2013) found sediment denitrification in restored oyster reefs in the Choptank River tributary of the Chesapeake Bay increased from April to June, peaked in August, and was lowest in November. Denitrification measurements were also taken at an oyster farm located several km nearer the river mouth 1 yr later, with peak denitrification in June, and a decline in August (Testa et al. 2015). There is a long-term water quality monitoring site (Chesapeake Bay Program Water Quality Data, Tidal Water Quality Monitoring Program station ET5.2; http:// data. chesapeakebay. net/ WaterQuality) between the 2 sites that records monthly chl a concentrations that we can use to roughly estimate peak chl a and determine if denitrification traces this pattern at these sites. For both locations, the highest denitrification measurements appear to follow the peak in annual surface chl a concentration. At the restored reef in 2010, both the chl a and denitrification peaks were in August, and at the oyster farm in 2011, peak denitrification was recorded 2 mo after the chl a peak in April. It appears that, similarly to Narragansett Bay, denitrification in these 2 oyster habitats tended to track water column productivity. The water column chl a concentrations in the Choptank River tributary were often more than twice as high as those in our study, which may explain why sediments were always net denitrifiers.
In Bogue Sound, 2 studies found that sediment denitrification in oyster habitats is lowest in winter and increases in spring to a summer peak before declining in fall (Piehler & Smyth 2011, Smyth et al. 2013b). However, a more recent study demonstrated this pattern does not hold across sites in the estuary (Onorevole et al. 2018). We could not locate any nearby water quality monitoring station data to com-pare with these studies, but chl a concentrations in nearby watersheds vary throughout the year and tend to be highest following peaks in monthly rainfall (Mallin et al. 1991). Taken together, it appears that across systems, rates of sediment denitrification in oyster habitats likely track water column productivity, if it is high enough to outpace oyster demand.

Ecological carrying capacity and oystermediated sediment nitrogen cycling
The long-term studies on Narragansett Bay's past and present N loading, alongside estimates of oyster populations at carrying capacity and past population sizes, present an opportunity to calculate the importance of oysters in the N cycle of Narragansett Bay under different oyster population and N loading scenarios. While the values we use here are specific to Narragansett Bay, the conclusions we make can be extrapolated to other coastal systems.
In Narragansett Bay, less than 1% of the historic oyster population remains (Zu Ermgassen et al. 2012), but there are efforts to reverse this trend through development of aquaculture and reef restoration. A recent model estimate suggests that oyster populations in Narragansett Bay could be 625 times larger than the 2010 population before the system reaches ecological carrying capacity (Byron et al. 2011). This estimate raises at least 2 important questions in the context of nutrient cycling: How will the Bay-wide N budget change if oyster populations increase 625 times? How did the N budget for Narragansett Bay look before the arrival of colonists and over-exploitation of oyster populations? Further, the influence of oyster biodeposition and oyster-driven changes in sediment denitrification and N 2 fixation may influence ecological carrying capacity and were not included in the model.
Using the values presented in Byron et al. (2011), Narragansett Bay at ecological carrying capacity would have approximately 5.93 g dry oyster tissue m −2 across the entire Bay (355 km 2 ). Assuming that oyster density is the same in the real world as in the carrying capacity model, we estimate that oysters currently influence 0.568 km 2 of sediment area in Narragansett Bay. From 2013 to 2015, the estimated annual N load to the Bay was 292 Mmol N (Narragansett Bay Estuary Program 2017), or 0.8 Mmol N d −1 if we assume uniform loading throughout the year. Using the sediment denitrification rates measured in this study and rates of denitrification in oyster digestive systems from a study conducted in nearby Ninigret Pond (Ray et al. 2019a), and assuming that these rates apply to any oyster habitat regardless of oyster density, we estimate that at current oyster populations, 671 mol N d −1 are removed via denitrification from oyster habitats in spring -approximately 0.1% of the N load (Table 3). At the modeled carrying capacity, this removal increases to 52% of the present daily N load or 142−175% of the daily prehistoric N load (assuming denitrification is constant despite lower N loading of 89−108 Mmol N yr −1 ; Nixon 1997).
In the summertime, we observed a near zero net N 2 -N flux but similarly high rates of NH 4 + recycling. Previously published estimates of summertime oyster metabolic N requirements (51.4 μmol N g DW −1 d −1 ; Langdon & Newell 1989) indicate that at carrying capacity, oysters in Narragansett Bay would require 108 kmol N d −1 -an amount equal to 37−44% of the daily prehistoric N load, or 14% of the modern N load. This high N demand (and seasonally important N removal via denitrification) requires the use of recycled N within the system, which oysters stimulate (Table 3). At carrying capacity, we estimate that oyster habitats recycle 985 082 mol N d −1 to the water column (Table 3), or 3.3−4.1 times the daily N load prior to colonization, and 1.2 times the modern N load, indicating that oysters drive efficient use of N in the system while simultaneously removing excess N. This recycled N can support additional phytoplankton growth, providing additional food for oysters.
Fall N 2 fixation in oyster habitats occurs at nearly the same rate as spring denitrification, with a current input of 605 mol N d −1 via N 2 fixation or 0.4 Mmol N d −1 at carrying capacity (Table 3). At carrying capacity, this N 2 fixation more than doubles the daily N load to prehistoric Narragansett Bay estimated by Nixon (1997), who suggested that N 2 fixation in the watersheds surrounding Narragansett Bay may have increased the prehistoric N load, but made no mention of N 2 fixation within the Bay itself. Fulweiler & Heiss (2014) did suggest that N 2 fixation may have been important in prehistoric Narragansett Bay. Nixon also suggested N from offshore may have been an important contributor to the high productivity of the Bay prior to anthropogenic effects. More recently, Oczkowski et al. (2016) estimated a prehistoric δ 15 N value for dissolved N in Narragansett Bay of 6 ‰ using samples from Native American shell middens. This value is close to the δ 15 N of dissolved N in the ocean, which is typically less than 7 ‰ (Chaves 2004). The δ 15 N associated with N from N 2 fixation is typically near 0 ‰, which, when combined with terrestrial N sources, could lead to a pooled δ 15 N of near 6 ‰. We propose that N 2 fixation in oyster habitats may have been an important seasonal contributor to the availability of N and productivity of prehistoric Narragansett Bay as well as other coastal ecosystems with large oyster populations and low nutrient inputs.
In systems with large oyster populations, keeping nutrients in the system is critical, and it is no surprise that there is typically little denitrification in oyster habitats in Narragansett Bay during periods of low productivity. Further, oysters may help to balance productivity across seasons by regulating sediment Denitrification can also occur in oyster digestive systems (Smyth et al. 2013a, Caffrey et al. 2016, Ray et al. 2019a). Here, we estimated N removal and recycling from oysters using published values from an oyster farm within 20 km of Narragansett Bay (Ray et al. 2019a) of 0.41 μmol N 2 -N ind. −1 h −1 and 1.42 μmol NH 4 + ind. −1 h −1 . Individual oysters in that study weighed 2.93 g DW, and we assumed rates would be the same across seasons Table 3. Daily Narragansett Bay-wide estimates of N fluxes from oyster habitats (mol N d −1 ) for current oyster populations ('current'), and at ecological carrying capacity ('capacity'). Removal values are calculated using net N 2 -N fluxes and recycling values using NH 4 + fluxes measured in this study. Negative removal values indicate net N fixation and addition of N to the system. We assumed the modern daily N load to equal 0.8 Mmol N d −1 (Narragansett Bay Estuary Program 2017) and the prehistoric daily N load to equal 0.27 Mmol N d −1 (Nixon 1997) denitrification and N 2 fixation. It is critical that model estimates of oyster carrying capacity consider and include the role oysters play in changing sediment N cycling. Byron et al. (2011) did not consider the impact of oysters on sediment N-cycling -either removal via denitrification or addition by N 2 fixationthough the role oysters play in regulating these processes is significant at an ecosystem scale. To be clear, we do not want to discredit Byron et al. (2011), as their estimates drove us to ask the questions presented in this manuscript. However, the evidence we present here clearly demonstrates that future models of coastal ecosystems with large bi valve populations should consider how these organisms regulate N-cycling.

CONCLUDING REMARKS
It is clear that sediment biogeochemical processes in oyster habitats are dynamic and may be dependent on water column productivity. We show that the net N 2 flux in oyster habitats switches between net removal via denitrification in the spring to net addition via N 2 fixation in the fall. This complexity may hinder accurate estimates of ecological carrying capacity and understanding of productivity and ecosystem structure in coastal habitats that once had significant oyster populations. Reducing N loads to coastal systems with large oyster populations may not necessarily hinder productivity due to high rates of N 2 fixation. As efforts to include oyster-mediated sediment denitrification in nutrient management planning are further developed, a holistic approach that considers weekly or even daily patterns of water column chl a and sediment OM quality may be necessary for accurate quantification.