Rebuilding A Collapsed Bivalve Population, Restoring Seagrass Meadows, and Eradicating Harmful Algal Blooms In A Temperate Lagoon Using Spawner Sanctuaries

Introduction

Bivalves are keystone organisms whose suspension feeding provides a wealth of ecosystem services including reducing phytoplankton biomass and other suspended particles (Officer et al., 1982; Hawkins et al., 1996; Barillé et al., 1997), controlling harmful algal blooms (Cerrato et al., 2004; Harke et al., 2011), cycling nutrients and organic matter between the benthos and the water column (Kautsky and Evans, 1987; Smaal and Prins, 1993), increasing light penetration (Newell and Koch, 2004), and facilitating the growth of submerged aquatic vegetation (Peterson and Heck, 2001; Carroll et al., 2008; Wall et al., 2008). In recent decades, as bivalve populations in many ecosystems around the world have declined due to overfishing, habitat loss, disease, and harmful algal blooms, these ecosystem services have been lost and there are currently few estuaries with natural densities of bivalves sufficient to exert ecosystem-wide effects (Newell, 2004; Lotze et al., 2006; Beck et al., 2011; Ermgassen et al., 2012). Such losses can have strong negative consequences for coastal ecosystems (Jackson et al., 2001; Lotze et al., 2006; Worm et al., 2006)

Reductions in bivalve filtration capacity leaves estuaries more vulnerable to environmental impacts associated with coastal human population expansion such as nutrient pollution (Lotze et al., 2006; Wall et al., 2011). Excessive nutrient loading accelerates primary production and can have a multitude of adverse impacts on affected ecosystems including the promotion of harmful algal blooms (HABs), hypoxia, acidification, and the subsequent loss of marine life and habitats (Valiela et al., 1992; de Jonge et al., 2002; Wallace and Gobler, 2021). There is evidence that these symptoms of eutrophication have been increasing (Breitburg et al., 2018; Anderson et al., 2021) and in many cases, these stressors can be exacerbated by climate change (Gobler et al., 2017; Breitburg et al., 2018; Griffith and Gobler, 2020).

During the past fifty years, the observed increase in eutrophication and decline in abundance of bivalve populations on Long Island and specifically Shinnecock Bay have been emblematic of global trajectories. During the 1970s, Long Island hosted the most prolific hard clam (Mercenaria mercenaria) fishery in the US (McHugh, 1983). Since that time and through the year 2012, landings of hard clams in the south shore lagoons of Long Island including Shinnecock Bay declined by more than 99% (NYSDEC, 1970-2011; Kraeuter et al., 2008) and densities of hard clams declined by more than 20-fold (Weiss et al., 2007; Kraeuter et al., 2008). In parallel with declining bivalve populations, Long Island has seen a 60% increase in the concentrations of nitrogen in groundwater discharging into estuaries (SCSWWP, 2020), a 90% loss of seagrass (NYSDEC, 2009), and a proliferation of HABs (Gobler et al., 2005; Gobler et al., 2008; Hattenrath et al., 2010; Hattenrath-Lehmann et al., 2013). Shinnecock Bay, in particular, has transformed from being an estuary that had never had a documented HAB prior to 1985 to experiencing paralytical shellfish poisoning blooms caused by Alexandrium catenella and associated shellfish bed closures and brown tides caused by Aureococcus anophagefferens on an annual basis (SCDHS, 1985-2021; Gobler and Sunda, 2012; Wang et al., 2012). As noted above, there were likely synergistic interactions between the loss of bivalves and the enhanced nutrient loads (Lotze et al., 2006) that made Shinnecock Bay particularly vulnerable to eutrophication and HABs, primarily due to its shallow nature (mean depth = 2 m; Cerrato et al., 2004).

Given the widespread declines of bivalve populations around the world (Jackson et al., 2001; Kraeuter et al., 2008; Beck et al., 2011), there is great interest in developing approaches to restore these populations (Duarte et al., 2020), including alleviating recruitment limitation via stocking substrate/shell material (oysters; Waldbusser et al., 2011), juveniles (oysters, clams; Peterson et al., 1995; Baggett et al., 2015), or adult spawning stock (scallops, hard clams; Kassner and Malouf, 1982; Doall et al., 2008; Tettelbach et al., 2011). Large quantities of small juvenile bivalves can be quickly reared and planted for restoration, but predation of bivalves is strongly size dependent, with rapid, mass mortality occurring for individuals smaller than predation refuge threshold sizes (Juanes, 1992; Peterson et al., 1995; Gosselin and Qian, 1997; Munroe et al., 2015). Consequently, restoration projects reliant on the distribution of small juvenile bivalves have met with limited success (Peterson et al., 1995; Geraldi et al., 2013). For free-spawning invertebrates, fertilization success is exponentially related to the density of individuals in a given area (Levitan, 1991; Mann and Evans, 2004) and conversely, overharvests and severe reductions in densities can lead to recruitment failure of bivalve populations (Peterson and Summerson, 1992; Levitan and Petersen, 1995; Peterson et al., 1996; Arnold et al., 1998; Liermann and Hilborn, 2001; Peterson, 2002; Kraeuter et al., 2005). Hence, an alternative approach to bivalve restoration has been enhancing spawning stock of reproductive adults to increase fertilization efficiencies and population reproductive output (Kassner and Malouf, 1982; Peterson et al., 1996; Peterson, 2002; Doall et al., 2008; Tettelbach et al., 2013; Tettelbach et al., 2015).In 2012 the Shinnecock Bay Restoration Program (ShiRP) was initiated to rebuild the hard clam populations of Shinnecock Bay, NY, USA. Prior to the start of this project, hard clam densities in Shinnecock Bay were at or below 0.7 hard clams m^-2 (Weiss et al., 2007) suggesting this population was recruitment limited (Peterson, 2002; Kraeuter et al., 2005; Bricelj, 2009). The primary strategy employed by ShiRP to restore hard clam populations was the establishment of spawner sanctuaries, or harvest free areas, planted with high densities of adult clams. These sanctuaries were intended to function as larval production zones, helping to increase the total reproductive output of hard clams in Shinnecock Bay through two primary mechanisms: 1) by increasing the total number of gametes spawned, and 2) by increasing fertilization efficiencies. The total number of gametes spawned was expected to increase with increases in the number of spawning stock added to the bay. Fertilization efficiencies were expected to increase with increases in localized clam densities, as the probability of successful fertilization is a function of the distance between males and females in free-spawning marine invertebrates (Levitan, 1991; Levitan and Petersen, 1995). Peterson (2002) suggested that the creation of spawner sanctuaries may be a particularly effective restoration strategy for M. mercenaria because of this species’ relatively long lifespan, and lack of reproductive senility (Peterson, 1983). Hard clam spawner sanctuaries may thus result in continued spawning and contributions to the larval pool for many years after transplantation, provided that the transplanted clams survive and meet the energy requirements necessary to ripen gonads and spawn. In addition, it is essential that the larvae be retained within areas of the estuary favorable to growth and survival, so that recruits can reach reproductive maturity and further promote population growth.In a sustained multi-year effort between 2012-2019, ShiRP planted over 3.2 million adult hard clams into spawner sanctuaries in Shinnecock Bay. In parallel, an intensive monitoring program evaluated the effectiveness of spawner sanctuaries at increasing hard clam recruitment and abundance and quantified the ecosystem-level impacts associated with changes in hard clam abundance and total biotic filtration. Monitoring of the density and reproductive condition of planted clams was used to determine whether clams were surviving and spawning, providing short-term measures of success each year. Analysis of data from annual benthic surveys of bay-wide hard clam recruitment and abundance, combined with analysis of commercial landings data from the New York State Department of Environmental Conservation (NYSDEC), provided medium term (3+ years) measurements of success as to whether hard clam populations were increasing in Shinnecock Bay. Bay-wide monitoring of hard clam larval abundance using digital-PCR, combined with hydrodynamic modeling, were employed to better understand larval temporal and spatial dynamics and transport within Shinnecock Bay, and determine if areas experiencing increased recruitment were receiving larvae produced by spawner sanctuaries. The monitoring of seagrass densities and general water quality (HABs, water clarity, chlorophyll a) across the estuary provided longer-term (6+ years) evaluation of broader effects of the restoration effort.

Methods

Locating Spawner Sanctuaries

Shinnecock Bay is located on the south shore of Long Island, NY (Figure 1A). It is part of a series of interconnected barrier-island lagoons that make up the South Shore Estuary Reserve, including, from west to east, Hempstead Bay, South Oyster Bay, Great South Bay, Moriches Bay, and Shinnecock Bay, which is the easternmost bay. Shinnecock Bay is microtidal (tidal range = 1 m), covers a ~40 km² area and has a constriction in the middle that separates the bay into distinct eastern and western portions (Figure 1B). An inlet to the Atlantic Ocean in the eastern half of the bay provides strong tidal flushing to this portion of the bay that attenuates to the west.

Figure 1 (A) Map depicting location of Shinnecock Bay, NY, USA, in the northeast US. (B) Close-up of the Shinnecock Bay denoting locations of the two spawning and harvest-free sanctuaries in red and the locations of sampling water quality and larval sampling stations (blue circles).

A series of criteria were set for the establishment of adult hard clam spawner sanctuaries in Shinnecock Bay. Regions that had previously been studied and shown to support robust clam growth, survival, and/or condition were clear priority locations. Regions that experienced extended periods of high temperatures (>25°C) and extreme levels of the harmful brown tide alga, Aureococcus anophagefferens (> 10⁶ cells mL^-1) were avoided due to adverse effects of these conditions on M. mercenaria (Ansell, 1968; Hamwii and Haskin, 1969; Greenfield et al., 2005; Gobler and Sunda, 2012). Sediments composed of extremely coarse sediments or anoxic, sulfidic muds, which can be harmful to bivalves (Wang and Chapman, 1999), was avoided. An exceedingly thick covering of macroalgae, which can be found in some regions of the Bay (Young and Gobler, 2016), was avoided assuring clams would not be smothered in seaweeds that might promote anoxia (Valiela et al., 1992). Regions with high levels of channel whelk (Busycotypus canalicultatus), a well-known predator on the south shore of Long Island (Polyakov et al., 2007), were also avoided. Given that hard clams have a two-week larval stage (Carriker, 2001), sites with extended residence times, located away from the Shinnecock Inlet were selected to ensure the retention of larvae within the estuary. Finally, sites selected based on these scientific criteria required regulatory and public approval.

Planting Hard Clam Spawner Sanctuaries

Within the spawner sanctuaries established in Weesuck Creek and Tiana Bay, adult clams were transplanted at high densities (mean ± SD = 27.1 ± 4.7 m^-2) onto numerous smaller plots (~1,850 m²) to increase fertilization efficiencies and the production of larvae, instead of being spread out at low densities across the entire area. In total, 30 GPS-delineated plots were established in Weesuck Creek, and 34 plots were established in Tiana Bay, spaced in a near equidistant manner from one another with the borders of each sanctuary area defined and GPS-delineated. Between 25,000 to 83,000 clams were planted in each plot yielding densities of 12.3 - 44.6 clams m^-2, with a mean planting amount and density of 50,300 clams and 27 clams m^-2, respectively. In total, 3,219,250 clams were planted from 2012 – 2019.

Sanctuary plots were planted with adult, market-sized (i.e., ≥ 25.4 mm shell width) clams purchased from seafood distributors. The clams were harvested from New York and Connecticut waters, primarily in the Long Island Sound where most of the current hard clam fishery exists for these states (Figure 1A). The clams were harvested by hand with clam rakes, sold to local distributors on the day of harvest, and stored under refrigeration by the distributors for up to three days prior to transport to vessels (R/V Peconic or R/V Shinnecock) for planting. To plant clams, the research vessel slowly (~ 1-2 knots) moved across the entire latitudinal and longitudinal confines of the sanctuary plot boundaries defined by the four corner GPS coordinates while clams were gently, evenly, and continually dispersed into the water from the deck of the pontoon vessels. Any visibly dead or broken individuals were removed to minimize the attraction of predators (Polyakov et al., 2007). Clams were planted during spring and fall when water temperatures were above 10°C to ensure clams would be active enough to dig into sediments. Planting ceased during summer months to avoid stressful temperatures (>25°C; Ansell, 1968; Hamwii and Haskin, 1969; Greenfield et al., 2005).

Monitoring Clams in Sanctuaries

Three plots in each of the two sanctuaries were surveyed annually during summer to quantify hard clam abundance and size. At each plot, divers hand excavated 20 randomly placed quadrats (0.25 m²) to a depth of approximately 15 cm. The shell height, length, and width of recovered clams were measured with calipers to the nearest 1 mm.

As a proxy for hard clam gametogenesis, clam gonad condition was monitored twice per year, near the beginning (June) and end (September) of the hard clam spawning season from 2013 - 2019. In 2020, gonad condition was monitored more frequently during summer (~weekly in June and July) to generate a refined understanding of the timing of spawning. On each sample date, ten clams were sampled from each of three plots in each sanctuary that were also monitored for survivorship (n=30 clams per sanctuary). The clams were either hand excavated by SCUBA divers (2012-2019) or collected using a clam rake from a boat (2020). Gonad ripeness (GR) was monitored qualitatively using a visual ranking system that ranged from 0 (no gonad) to 4 (most ripe), with declines in gonad ripeness providing an indication of spawning (Doall et al., 2008). Rankings were based on the relative amount of gonad tissue observed within the dissected visceral body, as well as the relative number and appearance of gametes observed in a tissue smear through a compound microscope as described by Doall et al., 2008.

Bay-Wide Benthic Surveys

Benthic surveys of ~60 randomly generated sites outside of the sanctuary plots were conducted by SCUBA annually during late summer from 2012-2021 to track changes in the abundance and size distribution of clam populations throughout eastern and western Shinnecock Bay. All sites were at least 200 m apart. At each site, divers excavated two quadrats (1 m²) to a depth of 15cm using a suction dredge to vacuum all quadrat contents into a 3-mm mesh collection bag. The contents of collection bags were sieved through a 3-mm screen on the boat, and all hard clams and crabs (mud crabs, Dyspanopeus sayi, blue crabs, Callinectes sapidus, lady crabs, Ovalipes ocellatus, spider crabs, Libinia emarginata, green crabs, Carcinus maenas, broad-clawed hermit crabs, Pagurus pollicaris, and long-clawed hermit crabs, Pagurus longicarpus) greater than this size were counted and measured for shell length, height, and width with calipers. Pre-market size clams (<25.4 mm shell width) were binned into one of three size categories (shell width) based on size-age relationships established for clams in the hydrologically-connected lagoon, Great South Bay (Kraeuter et al., 2005): 1) < 9 mm (0-year-olds); 9-17 mm (1-year-olds); 18-25 mm (2-year-olds). Market size clams were binned into four size categories corresponding to the commercial size classes of littleneck (26-35 mm), topneck (36-38 mm), cherrystone (39-41 mm), and chowder (>41 mm). The size ranges for each commercial size category were adapted from Kraeuter et al. (2005), with the addition of the ‘topneck’ size category added between littlenecks and cherrystones to correspond with NYSDEC commercial landings data.

Quantification of Hard Clam Larvae

Digital PCR was used to detect and quantify the concentration of hard clam larvae in water samples collected via twice weekly cruises across Shinnecock Bay during June and July 2020, when hard clams are known to spawn on the south shore of Long Island (Doall et al., 2008). A total of 14 cruises were performed, two per week, from 6/10/2020 – 7/27/2020, sampling from seven stations spread across eastern and western Shinnecock Bay, including: Far East (40.86386, -72.44838), East Central (40.8659, -72.4718), Bridge East (40.85854, -72.49019), Bridge West (40.84375, -72.51406), Tiana (40.85523, -72.5381), Weesuck (40.83698, -72.5678), and Far West (40.81905, -72.58297) (Figure 1). At each site, two replicate 20L water samples were obtained from ~0.5m below the surface and were concentrated onto a 20µm sieve, the contents of which were backwashed into two, 50mL centrifuge tubes vials that were immediately placed on ice and stored at -20°C within two hours of collection.

To prepare samples for digital PCR analysis, total cellular genomic DNA was extracted using the Qiagen Blood and Tissue Kit (Qiagen, Germantown, MD, USA). Samples were filtered onto 0.2 µm polycarbonate filters that were transferred to a 5 mL bead tube and treated with a lysis buffer and proteinase K to chemically lyse all cells and remove non-DNA organic contents. The lysate was transferred to a new tube and 50% v/v solution of ethanol was added to precipitate total genomic DNA. The mixture was subsequently passed through a mini-spin column via centrifugation (13000 g) to capture the DNA. The columns were further washed with a buffer containing ethanol, followed by elution of high-quality DNA with 100 µl nuclease-free water. The eluted samples were analyzed on a Nanodrop Spectrophotometer (Thermo Scientific, Waltham, MA, USA) to ensure nucleic acid recovery and quality. The purified DNA samples were stored at -80°C until digital polymerase chain reaction (dPCR) analysis.

Digital PCR analysis was conducted using the chip-based Applied Biosystems™ QuantStudio™ 3D Digital PCR System (Applied Biosystems, Foster City, CA, USA) to quantitatively identify the number of mitochondrial cytochrome c oxidase subunit I (COI) gene copies from M. mercenaria. To identify unique regions of the COI gene for this species, sequences for various bivalve species found across the US East Coast, including the closest relative Gemma gemma, were obtained from publicly available databases, NCBI GenBank and Barcode of Life Data Systems (BOLD; Ratnasingham and Hebert, 2007). The sequences were aligned, and unique regions were used to design forward (GGGCTCTTTCTCATTCTGGTAG), and reverse (AGCAACACCATAACACCCG) primers, and probe (TGGGTGGTGCATCTTCTATTTTGGCG). Each assay was validated and optimized using the dPCR system prior to sample analysis using synthetic double-stranded DNA fragments of the target genes as standards (gBlocks, Integrated DNA Technologies, Coralville, IA, USA). Lyophilized gBlocks were resuspended in 25 μL of IDTE buffer + 100 ng/μL polyA carrier (Roche Life Sciences, South San Francisco, CA, USA) used to increase the recovery of synthetic standards (Miyaoka et al., 2016), quantified using a Qubit 4 Fluorometer (Thermo-Fisher, Waltham, MA, USA), and serially diluted to prepare standards with final concentrations of 800 copies µL^-1. Optimization trials, testing gradients of annealing temperature, primer-probe concentrations and numbers of cycles were conducted to identify optimal thermocycling conditions for the assay.Digital PCR amplifications were performed in 14.5 µL reaction mixtures consisting of 7.25 µL of QuantaStudio 3D digital PCR Master mix v2 (2x stock solution), 0.725 µL Taq Man assay primer and probe mix (20x stock solution, see Table 5.1 for final concentrations), 1.525 µL nuclease free water and 5 µL sample DNA. All samples were originally analyzed using a maximum 5 µL of extracted DNA to achieve an on-chip concentration in the optimal range of 200-2000 c/µL. The dPCR reactions were loaded onto QuantStudio™ 3D Digital PCR Chip V2 chips containing 20,000 well partitioning with the QuantStudio™ 3D Digital PCR Chip loader (Applied Biosystems, Foster City, CA, USA), sealed with immersion fluid and the chip lid per the manufacturer’s instructions. All chip preparation was performed in less than one hour per manufacturer’s recommendations to minimize degradation. Loaded chips were amplified using a ProFlex™ 2x Flat PCR System thermocycler (Applied Biosystems, Foster City, CA, USA). Thermocycling conditions included initial denaturation at 95°C for 10mins, followed by 40 cycles of three steps: 95°C for 10mins, 95°C for 10mins and 95°C for 10mins. Amplified chips were brought to room temperature to prevent condensation before imaging on the QuantStudio™ 3D Digital PCR instrument (Applied Biosystems, Foster City, CA, USA). All samples were run as duplicates, along with a negative (nuclease free water) and positive (gBlock standards, 800 copies µL^-1 concentration) control.Imaging data derived from the QuantStudio™ 3D Digital PCR instrument were analyzed using the Applied Biosystems QuantStudio^® 3D AnalysisSuite™ cloud software. This software provided quality control steps on a per chip basis determining wells suitable for further analysis. In this study the default quality threshold of 0.5 was used for all chips. Chips were also manually inspected for equal distribution of positive wells across the chips and chip damage, such as large bubbles or evaporation, resulting in loss of readable wells in which chips were omitted and the sample was reanalyzed. Software derived fluorescence (call) thresholds delineating the unamplified wells (negative calls) and amplified wells (positive calls) were manually reviewed for each chip and adjusted to a common threshold per assay based on the ranges of the positive control and negative control clusters. Additionally, the spread of reads along the secondary assay (non-target dye) was manually reviewed in which wells identified as positive located largely outside the range of the positive control clusters on the secondary axis were identified as no amplification to reduce false positives. The negative and positive well counts were then converted to absolute quantification (copies µL^-1) by the software using Poisson statistics and corrected for dilution/concentration factors during sample collection (filtration), DNA extraction, and PCR reaction preparation. Sample concentrations are reported in M. mercenaria COI gene copies per liter.

Water Quality Sampling

Discrete water samples were collected weekly-to-monthly from May through October from selected locations in Shinnecock Bay from 2012-2021, with weekly sampling in summer (June – September) and monthly in spring and fall. At each station, temperature and salinity were measured at the surface and near the bottom using a YSI^® 556 sonde to confirm that the water column, which was typically 2 m, was well-mixed (Wilson et al., 1991). Water samples were collected in 20-L carboys at ~0.5 m at each station and returned to the laboratory for further analyses within 30 minutes. Triplicate chlorophyll a samples were collected on GF/F glass fiber filters (pore size = 0.7 µm), frozen, and analyzed by standard fluorometric methods (Parsons et al., 1984). Densities of the brown tide alga, Aureococcus anophagefferens, were quantified on a flow cytometer using a species-specific immuno-assay on water samples preserved in 1% glutaraldehyde (Stauffer et al., 2008).

A 20-year (2000-2020) dataset compiled by the Suffolk County (NY, USA) Department of Health Services (SCDHS) water quality monitoring program was analyzed to compare conditions before and after the creation of hard clam spawner sanctuaries. This program samples dozens of sites across the south shore of Long Island on a biweekly-to-monthly basis. Sampling occurred monthly and during the day without regard for tidal state. Sites are evenly distributed across Shinnecock Bay, typically separated by 1 – 3 km (Figure S1). Chlorophyll a and nutrients were measured using standard oceanographic methods (Parsons et al., 1984), and Secchi disk depths were quantified. Water quality conditions were compared pre- and post-hard clam spawner sanctuaries and over time.

Seagrass Assessment

New York has lost 90% of its seagrass meadows since the 20^th century (NYSDEC, 2009). Here, classification of orthoimagery was performed to generate an estimate of eelgrass (Zostera marina) coverage change from 2013 to 2017. The imagery was obtained from the USDA’s National Agriculture Imagery Program (NAIP) which has collected high-resolution four-banded aerial images of New York state every other year since 2011. NAIP imagery has a 1-meter ground sample distance resolution, meaning a raster cell size of 1m by 1m. The NAIP imagery is collected to survey agricultural land during peak growing conditions (June – August) which fortuitously aligns with the peak eelgrass growing season on Long Island, NY. Since the intention of the imagery collection does not include surveying coastal waters, poor sea state conditions and sun angles can obscure the benthos. The 2011 and 2019 imagery were completely obscured and unusable, the 2021 imagery is not yet available, and therefore, the 2013 and 2017 orthoimagery were used for this analysis.

The western side of Shinnecock Bay, where spawner sanctuaries were established, was selected as the target region to assess eelgrass coverage change. Classification of eelgrass coverage change was partly constrained by the 2013 imagery as sun glare obscured the imagery on the western and easternmost edges of the bay, further reducing the viable area for analysis (Figure S2). Depths beyond two meters were excluded as target regions as 99% of eelgrass in NYS is shallower than this depth (NYSDEC, 2009) and deeper regions become unusable for imaging of benthic structures.Natural color bands from the NAIP imagery (bands 1-3, red, green, blue) were used to classify bay bottom into four categories (shallow non-seagrass substrate, deep non-seagrass substrate, shallow seagrass, and deep seagrass) through a supervised object-based support vector machine via Esri’s ArcGIS Pro 2.9.1 Image Classification Wizard tool following specifications for object creation described by Davenport et al. (2017). Forty training samples were chosen for each class with a minimal pixel coverage of 20px by 20px. After model training and classification, the deep and shallow seagrass classes were combined, extracted, and converted to a polygon feature class to allow for geometry calculations. A union was performed on the 2013 and 2017 seagrass polygons. Areas of overlap in eelgrass coverage between the two years were assigned as no change for eelgrass coverage, areas with eelgrass only from the 2017 polygon were assigned as a gain in eelgrass coverage, and areas from solely the 2013 polygon were assigned as losses in eelgrass coverage. The total area was then calculated for each of the three assigned categories of eelgrass coverage. Ground-truthing was performed via comparison of model outcomes to diver surveys to affirm Z. marina coverage in the predicted locations.

Bio-Optical Model

A previously developed bio-optical model based on the absorption and scattering of light by specific components of coastal waters (Gallegos and Kenworthy, 1996; Gallegos and Biber, 2004 2008) was calibrated for Shinnecock Bay and the characteristics of the dominant seagrass in Shinnecock Bay, Zostera marina. The optical model is based on the properties of three optically active constituents: non-algal particulates (TSS), phytoplankton (chlorophyll a), and colored dissolved organic matter (CDOM). Chlorophyll a was quantified as described above, TSS was determined gravimetrically (Parsons, 2013) and CDOM was measured by filtering seawater through a 0.7µm glass fiber filter and measuring levels on a Trilogy Laboratory Fluorometer with a CDOM module with ultraviolet GUI selection. Thirteen sites across Shinnecock Bay (Figure S3) were sampled during the day at peak water temperatures (August), a time that may be stressful for the seagrass, Z. marina (Hammer et al., 2018). Measurements made in August 2017 were compared to August 2021. Previous studies have demonstrated that light levels below a minimum physiological requirement (typically 15 – 40% of incident surface light) result in depth-limits to seagrass (Kenworthy and Fonseca, 1996; Steward et al., 2005). Therefore, the calibrated optical model was inverted to set threshold concentrations for water quality parameters to meet the water column light requirement for Z. marina at a given water depth (Gallegos and Moore, 2000).

Hydrodynamic Model

The hydrodynamic model, Environmental Fluid Dynamics Code (EFDC; Hamrick, 1992) was constructed and optimized for Shinnecock Bay to estimate the movement of water and bivalve larvae in this system. The model solves the three-dimensional, time-dependent equations of motion for variable-density water. An equation of state, which determines water density as a function of temperature and salinity, is used to couple conservation of mass (salinity) and energy (temperature) equations to the hydrodynamic equations of motion. Vertical mixing due to turbulence is simulated using the Mellor-Yamada turbulence closure scheme (Mellor and Yamada, 1982) as modified by Galperin et al. (1988). Time-dependent boundary conditions in EFDC can be specified at inflow and outflow boundaries. A complete description of the EFDC model, including numerical algorithms is provided by Hamrick (1992) and Hamrick (1996).

For this study, the spatial extent of the model included Shinnecock Bay, the Atlantic Ocean, and Moriches Bay (Figure S4). Inputs to the hydrodynamic model were obtained from NOAA tidal gauges in Moriches and Shinnecock Bays at their connections to the Atlantic Ocean and weather data from the Francis S. Grabreski Airport. Bathymetric data were obtained from the NOAA Coastal Relief Model (NOAA, 1999) and sounding surveys. The hydrodynamic model grid comprises 901 total cells. The aerial extent of the grid is about 40 km (25 miles) in the east-west direction, and 5 km (3 miles) in the north-south direction (Figure S4). The average dimensions of the cells are 370m (east-west) and 200m (north-south). The bays are relatively shallow (1 – 4m), and, in the absence of any observed vertical stratification, were simulated as a single vertical layer. The hydrodynamic model was calibrated to water surface elevations and temperature measured by NOAA at the Ponquogue Point station in Shinnecock Bay and the Moriches CGS station in Moriches Bay. Temperature data were obtained from the Suffolk County monitoring stations in Shinnecock, Quantuck, and Moriches Bay (SCDHS, 1985-2021).In the first iteration of the model, a conservative tracer (i.e., dye) was released from each spawner sanctuary to visualize water dilution and movement on the time scale of the larval life cycle (two weeks; Carriker, 2001) using physical conditions (winds, tides, temperature, salinity) found for Shinnecock Bay during summer when M. mercenaria spawns. In a more refined iteration, the transport of neutrally buoyant particles that mimic clam larvae was calculated using velocity vectors based on a grid resolution of approximately 250 meters and a time interval of 5 minutes. Discrete particle locations were computed independently of the hydrodynamic model grid extents and used the velocity magnitude and direction from the nearest model cell. This type of analysis provides hypothetical particle tracks that may vary depending on grid size, computation time interval, release locations, release times, and complexity of circulation patterns. Particle movement was tracked using the tidal and wind conditions present during July 2020, when a high abundance of larvae was quantified in Shinnecock Bay. Particle movement was determined for both spawner sanctuaries and for eastern Shinnecock Bay, with one particle placed in each grid cell for each region (14 in Weesuck Creek, 17 in Tiana Bay, and 3 in eastern, central Shinnecock Bay).

Hard Clam Landings

The commercial landings of hard clams from Shinnecock Bay and near-by lagoonal estuaries were tracked using long term data sets publicly available from the New York State Department of Environmental Conservation (NYSDEC). In New York State, the legal harvest size limit for hard clams is one inch (25.4 mm) shell width, which represents three years of growth (Kraeuter et al., 2005). Landings are reported to the NYSDEC by shellfish distributors, including information on quantity, date of harvest, harvest area, and clam size. Clam size is classified into one of four commercial size categories, including, from smallest to largest, littleneck, topneck, cherrystone, and chowder, with cherrystone and chowder size clams often reported together as a “large mix”. Prior to 2012, the NYSDEC tracked only total hard clam landings of all size categories combined, with landings expressed in terms of bushels. From 2012-2020, more detailed landings data is available from the NYSDEC, including information on the number of individuals harvested in each size category. To convert between individuals and bushels, the NYSDEC uses the following conversion factors: one bushel equals 600 littlenecks, or 200 topnecks, or 150 cherrystones, or 100 chowders, or 120 large mix (cherrystone and chowder mixed). Harvest locations are reported according to specific areas delineated by the NYSDEC, and from 2008 to 2021 separate harvest areas were delineated for western (SS12) and eastern (SS13) Shinnecock Bay. Prior to 2008, landings from the entire Shinnecock Bay and the adjoining smaller lagoon Quantuck Bay were grouped into one harvest area.

Ecosystem Filtration Rates

To quantify the potential filtering effect of M. mercenaria populations in Shinnecock Bay, biological filtration rates of clams within eastern and western Shinnecock Bay were determined from 2014 - 2021. Filtration capacities were derived from allometric relationships established between M. mercenaria weight and clearance rates for larger clams (>16 mm) fed natural seston (Hibbert, 1977; Doering and Oviatt, 1986) and for smaller clams (<16 mm) fed natural seston (Grizzle et al., 2001). Application of these relationships to previously reported clearance rates of M. mercenaria feeding on natural seston from Shinnecock Bay (Wall et al., 2008) yielded rates consistent with those measured empirically in that study. Size-specific clearance rates from the allometric relationship were applied to the densities of each surveyed size class of clam surveyed in each side of the bay and the clams planted in sanctuaries each year. These basin-wide clearance rates were further applied to the total water volumes of eastern and western Shinnecock Bay as determined via GIS and the NOAA Coastal Relief Model (NOAA, 1999) to generate basin-wide turnover times.

Results

Establishing Sanctuary Areas

Two areas within western Shinnecock Bay were chosen for the location of hard clam spawner sanctuaries: Tiana Bay and Weesuck Creek (Figure 1B). Prior studies documented high survivorship and condition index of clams planted in these regions (Padilla and Doall, 2012). Hydrodynamic models indicated that dye released from Weesuck Creek and Tiana Bay in western Shinnecock Bay dispersed to regions east and west of these sites over the two-week larval period of M. mercenaria (Figure 2). While both sites had a large percentage of the dye retained in western Shinnecock Bay, dye from Tiana Bay mixed across nearly all of eastern Shinnecock Bay and dye released within Weesuck Creek mixed across a portion of eastern Shinnecock Bay (Figure 2), suggesting larvae spawned in these sites would be spread across, and retained within, the Bay. SCUBA surveys found sanctuary regions to be comprised of appropriate sediment types (sandy mud), with moderate levels of macroalgae (<25% coverage), and an absence of channel whelk. Water quality data (SCDHS, 1985-2021) revealed these locations experienced more moderate temperatures (<25°C) and harmful brown tides (<10⁶ cells mL^-1) than western regions in this estuary. Following approval by the NYSDEC, Southampton Town Trustees and Southampton Bayman’s Association, precise areas within Tiana Bay and Weesuck Creek were designated as harvest-free sanctuaries in 2012 (Figure 1B).

Figure 2 Portion (0 – 1.0) of modeled dye remaining 14 days after release from (A) Weesuck Creek and (B) Tiana Bay. Arrows depict tidal flow speeds at the end of 14 days.

Clam Plantings

Adult hard clam plantings began in 2012. By the spring of 2015, more than one million clams had been planted (Figure S5A). Two and three million clams were planted by 2016 and 2017, respectively, while a slower addition of tens to hundreds of thousands of clams were planted annually from 2018-2019 (Figure S5A). By the end of 2019, 1,540,200 clams had been planted into 30 plots in the Weesuck Creek spawner sanctuary, and 1,679,050 clams had been planted into 34 plots in the Tiana Bay sanctuary, at densities ranging from to 12.3 to 44.6 clams m^-2 (Figures S5B,C). Clams planted were littleneck size and larger (>25mm).

Survivorship and Spawning of Planted Clams

In Weesuck Creek and Tiana Bay, hard clam densities in spawner sanctuaries have remained consistent across years (Figure 3A). Almost all hard clams collected during surveys were large adults (>39 mm shell width), with low juvenile abundance indicating low recruitment within the sanctuary areas (Figure S6). Therefore, the steady hard clam population size at the planting sites indicated high survivorship of the adults that were planted, rather than new recruitment replacing losses due to mortality. In addition, a decrease in mean gonadal ranks during the hard clam spawning season (summer) indicated annual spawning activity within spawner sanctuaries (Figure 3B). Combining years and sites, mean gonad ranks in late spring (2.85 ± 0.28; mean ± SD) were significantly higher than in late summer (2.24 ± 0.28) (t=4.93; p<0.001; paired two sample t-test for means; Figure 3B), with summer declines in gonadal rank averaging 21%. To further resolve the timing of spawning during summer, weekly monitoring was performed in 2020 which revealed high GR from late May through early July, but a large decline in GR thereafter, from 3.9 in early July to 2.7 by mid-July in Weesuck Creek, and from 3.7 to 2.9 in Tiana Bay, changes indicative of spawning (Figure 4A).

Figure 3 (A) Density of adult hard clams in spawner sanctuaries in Weesuck Creek and Tiana Bay over time shown as box plots where the whiskers are maximum and minimum values, the center line is the median and the boxes encompass the upper and lower quartiles. (B) Mean ( ± SE) gonad rank of hard clams in spawner sanctuary plots (n=3) in Weesuck Creek and Tiana Bay at the start and end of the spawning season, from 2013 to 2020, with blue shaded regions depicting the spawning season and the tan shaded region depicting the period of gametogenesis.

Figure 4 (A) Mean ( ± SE) gonad rank of clams in spawner sanctuary plots (n=3) during the summer of 2020 and (B) Concentrations of hard clam larvae in the water column (B) in Shinnecock Bay during summer 2020 (6/10/20-7/29/20). See Figure 1 for details of locations.

Spatial-Temporal Distributions of Hard Clam Larvae

Quantification of hard clam larvae across seven sites in Shinnecock Bay revealed low levels (< 1 x 10⁴ gene copies L^-1) through June and early July, with small peaks within the eastern bay on June 10 and at isolated locations in early July (Figure 4B). On July 17^th, the highest density of hard clam larvae to that date (5.8 x 10⁴ copies L^-1) was quantified within Tiana Bay with a smaller peak level found at the Far East station (2.3 x 10⁴ copies L^-1; Figure 4B). The timing of this peak corresponded with the timing of the largest decrease in gonadal rank of clams in the spawner sanctuaries (Figure 4A). Three days later (July 20), densities similar to those found in Tiana Bay on July 17^th were quantified to the east at the sites east and west of the Ponquogue Bridge in the center of the bay (6.2 x 10⁴ and 5.0 x 10⁴ copies L^-1, respectively), with high levels remaining in Tiana Bay but also appearing in the center of Eastern Shinnecock Bay (East Central; ~3 x 10⁴ copies L^-1 for both locales; Figure 4B). Visual examination of formalin-fixed plankton samples from net tows confirmed the presence of high concentrations of bivalve larvae at the Ponquogue Bridge on July 20. Hard clam larvae concentrations declined, but remained at moderate levels, through late July with 1.2 x 10⁴ copies L^-1 at the East Central site on July 24 and an increase in the Far Eastern site on the same date (Figure 4B).

To mimic the sudden increase and presumed release of larvae from the spawner sanctuaries on July 17^th, the hydrodynamic model was used to release neutrally buoyant particles based on tidal and wind conditions present on July 17, 2020. Hydrodynamic transport conditions near the mouth of spawner sanctuaries were relatively small in a north-south direction moving from the mouth out into Shinnecock Bay where velocities were relatively higher and predominantly in the east-west direction with the highest velocities occurring along the southern side of Shinnecock Bay. Over a two-week period, particles migrated east from Tiana Bay and Weesuck Creek, past the Ponquogue Bridge, and settled in eastern Shinnecock Bay, with the majority of Tiana Bay-released particles migrating to eastern Shinnecock Bay and Weesuck Creek-released particles displaying both eastern and western dispersal (Figures 5A, B). In contrast, particles released from eastern Shinnecock Bay were tidally transported into the Atlantic Ocean in four days (Figure 5C).

Figure 5 Hydrodynamically modeled 14-day dispersal of particles released from (A) Weesuck Creek, (B) Tiana Bay, and (C) eastern Shinnecock Bay. Yellow stars indicate starting points of particles released. Model runs used tides and wind directions and speeds starting on July 17^th, 2020. Release points matched the number of model grid points at each sanctuary (Supplementary Figure 2). Dark blue points are particle trajectories in the first week whereas light blue points are locations in the second week. All particles from eastern Shinnecock Bay were irreversibly advected to the ocean after two days and thus particles tracks are shown for those days only.

Surveys of Hard Clam Populations, 2012-2021

Benthic surveys of adult hard clams across the non-sanctuary zones of Shinnecock Bay revealed increases in densities from 2012 – 2021 with distinct patterns for the eastern and western Bay (Figure 6A). When benthic surveys began in 2012 in western Shinnecock Bay and 2015 in eastern Shinnecock Bay, clam densities were 1.2 and 0.4 clams m^-2, respectively, with populations being comprised primarily of large individuals > 36 mm (Figure 6A). In the subsequent years, both locations experienced steady increases in densities caused almost exclusively by smaller size classes of hard clams (Figure 6A). In the western half of the bay, clam density increased 50% between 2012 and 2021, with average densities of 1.8 and 1.9 clams m^-2 found in 2020 and 2021, respectively. In the eastern half of the bay, increases in clam densities were more dramatic, with increases occurring from 0.4 clams m^-2 in 2015 to 1.3 clams m^-2 in 2016, to 1.7 clams m^-2 in 2018, to 3.0 clams m^-2 in 2019, 3.5 clams m^-2 in 2020, and 7.3 clams m^-2 in 2021. In total, clam densities increased 18-fold from 2015 to 2021 in eastern Shinnecock Bay with a clustering of sites across the center of the eastern Shinnecock Bay displaying 6 – 36 clams m^-2 by 2021 (Figure 6B). The higher recruitment in eastern Shinnecock Bay has resulted in a dramatic shift in the population size structure (Figure 6A). For example, large clams (>36 mm; topnecks, cherrystones, and chowders) represented 60% of the population in the eastern Bay in 2015 but were <10% of the population by 2021 as small clams (<35 mm; littlenecks and juveniles) increased from 40% to 90% of the population (Figure 6A). The most dramatic change was among smallest size classes of clams (< 17mm) that were not detected in 2015 but comprised more than a one quarter of clams surveyed in 2020 and 80% in 2021 (Figure 6A). In western Shinnecock Bay, a smaller shift in population size structure was observed, with the dominance of small juvenile clams (<17 mm) increasing from 2% in 2012 to 25% in 2021.

Figure 6 (A) Density of hard clams in eastern and western Shinnecock Bay over time. Size classes are divided by shell width and include 0 years old (< 9 mm), 1 year old (>9 - < 16), 2 years old (>17 - < 25), littlenecks (> 25 - < 35 mm), topnecks (< 35 - >38 mm), and cherrystones + chowders (>38 mm). (B) Kriging interpolation of hard clam densities in Shinnecock Bay in 2021.

Commercial Landings of Hard Clams

Commercial landings provide another indicator of recent hard clam population changes in Shinnecock Bay. Annual hard clam landings from Shinnecock Bay have increased since the onset of the planting of adult clams in 2012 by more than sixteen-fold from 878 bushels in 2012 to 13,983 bushels in 2021 (p<0.0001; linear regression with time), with a modest rate of increase from 2012 – 2015 and a more rapid rise from 2016 – 2021 (Figure 7A). This increase has been stronger in the eastern half of Shinnecock Bay than the west (Figure 7A) and collectively the landings in 2021 were at a level (>14,000 bushels) not seen in Shinnecock Bay and surrounding estuaries since 1981 (Figure 7B). This sharp rise in landings has been regionally specific to Shinnecock Bay as landings in other, interconnected south shore Long Island lagoons have declined or had much smaller increases during the same time period (Figure S7). In the neighboring Moriches Bay, for instance, landings declined by 42% from 2012 to 2021 (Figure S7). While landings in Great South Bay increased from 2015 – 2021, landings in Shinnecock Bay now exceed those in Great South Bay, which is approximately ten-times the aerial size of Shinnecock Bay (Figure 1; Figure S7).

Figure 7 Hard clam landings in Shinnecock Bay as reported by the NYSDEC. (A) Landings in eastern and western Shinnecock Bay since 2008, (B) Landings in all of Shinnecock Bay and the adjoining smaller lagoon Quantuck Bay since 1970, (C) Landings in all of Shinnecock Bay since 2012 for three sizes of hard clams displayed as proportions, (D) Landings in all of Shinnecock Bay since 2012 for three sizes of hard clams.

Among the increased landings within Shinnecock Bay, the increase has been disproportionately among smaller clams. From 2012 to 2020, annual landings of littleneck clams increased faster (26-fold increase) than annual landings of larger cherry and chowder clams (seven-fold increase), indicating that the total increase in landings was largely due to recent population increases rather than solely increased fishing pressure (Figure 7C). With the disproportionate increase in landings of littleneck clams, the percentage of littleneck clams in total landings increased from 45% in 2012 to 67% in 2020 (Figures 7C, D).

Total Hard Clam Population and Ecosystem-Based Clearance Rates

The combined clams planted, and clam densities measured across the Shinnecock Bay provides an estimate of the total number of hard clams in eastern and western Shinnecock Bay (Figure 8A). In western Shinnecock Bay, this number rose significantly with time (p<0.01) from 1.9 x 10⁷ clams in 2014 to 4.2 x 10⁷ clams in 2021 whereas in the eastern portion of the bay the population rose from 0.6 x 10⁷ clams in 2015 (first year of full survey) to 1.1 x 10⁸ clams in 2021 (p<0.01; Figure 8A). Increased population sizes on each side of the bay translated into a more rapid ecosystem clearance rate and thus shorter biological turnover time of the bay over time when applying allometrically determined, size-specific clearance rates of hard clams to the in-situ population (Figure 8B). The western portion of Shinnecock Bay had a turnover time of three weeks in 2014 but this rate dropped to 11 days by 2017 and dropped to 10 days in 2020 (Figure 8B). Changes in the deeper, eastern portion of the bay were more dramatic with a biological turnover time of 104 days in 2015 that progressively dropped to 18 days by 2020 (Figure 8B). Applying similar calculations to the sub-bay regions hosting hard clam spawner sanctuaries indicates biological turnover rates in Weesuck Creek and Tiana Bay were initially two-to-three weeks and are now both ~ three days.

Figure 8 (A) Total number of hard clams in Shinnecock Bay since 2014 (landings plus sanctuary plantings), (B) Turnover times of eastern and western Shinnecock Bay based on quantified hard clam densities and allometrically estimated clearance rates.

Changes in Water Quality

Since the commencement of hard clam restoration in Shinnecock Bay, there have been distinct changes in several water quality parameters. Brown tides caused by A. anophagefferens had been annual events on the south shore of Long Island including Shinnecock Bay from 1985-2016 (Supplementary Table 1). From 2007 – 2016, brown tides with peak densities from 0.3 – 1.9 x 10⁶ cells mL^-1 had occurred annually in this system with particularly dense and sustained blooms detected in 2012 and 2016 (> 1.5 x 10⁶ cells mL^-1; more than six weeks; Figure 9). Since 2016, a brief, small brown tide bloom (~ 1 x 10⁵ cells mL^-1; two weeks; Figure 9) occurred in 2017, with no brown tides from 2018 – 2021 (<10⁴ cells mL^-1; Figure 9). Since 2012, there has been a significant decrease in the densities of A. anophagefferens and chlorophyll a concentrations measured in western Shinnecock Bay (stations 180 and 190; p=0.002 and p=0.0001; linear regressions with time) during the months of April through November, when clams are actively feeding (Figures S8A,B). There has also been a significant increase in the Secchi disc depth of western Shinnecock Bay (p=0.03; linear regression with time; Figure S8C). In contrast to Shinnecock Bay, other interconnected south shore lagoons not experiencing restoration have experienced significant increases in chlorophyll a and harmful brown tides caused by A. anophagefferens, as well as a significant decrease in water clarity during the same time period (Gobler et al., 2019).

Figure 9 Maximal densities of brown tide (Aureococcus anophagefferens) in Penniman Creek and Weesuck Creek, western Shinnecock Bay, from 2012 – 2021.

Light and Seagrass

Digital aerial imagery of seagrass in the western portion of Shinnecock Bay revealed that, among the areas visible via aerial data, there was an expansion of aerial coverage from 2013 through 2017 of 5.75 x 10⁵ m², as well as loss of 1.01 x 10⁵ m², yielding a net gain of 4.74 x 10⁵ m² (Figure 10A). Gains in seagrass were primarily concentrated in the western extent of the surveyed region, but also occurred within and near the Tiana Bay spawner sanctuary (Figure 10A). Sampling of parameters used in the bio-optical model demonstrated that there has been a substantial increase in light availability in western Shinnecock Bay between 2017 and 2021 that has been driven by the reduction in chlorophyll a in this portion of the bay (Figure 10B; Table S2). Light levels within the far western region and specifically in the Weesuck Creek and Tiana Bay spawner sanctuaries have been estimated to increase by up to 10% while light levels in some regions to the east have decreased, potentially associated with the transport, deposition, and resuspension of sediments recently deployed to fortify the barrier island near the ocean inlet.

Figure 10 (A) Change in Z. marina coverage in western Shinnecock Bay, 2013-2017; (B) Bio-optically modeled changes in light levels in Shinnecock Bay, 2017-2021.

Discussion

The twentieth century has witnessed the collapse of suspension feeding bivalves populations across the globe often with severe negative ecosystem consequences (Jackson et al., 2001; Cerrato et al., 2004; Lotze et al., 2006; Beck et al., 2011). In response, there have been a multitude of efforts to restore these populations using a variety of approaches (Duarte et al., 2020). This study demonstrated that planting dense aggregations of adult hard clams (M. mercenaria; 27 clams m^-2) in regions supportive of high clam survival, annual spawning, and maximal larval retention led to a progressive reestablishment of the population. Since 2012, the start of restoration activities, densities of hard clams in western and eastern Shinnecock Bay rose by 50% and 1,700%, respectively. Concurrently, landings of hard clams rose by 16-fold, biological turnover times decreased, levels of harmful algal blooms and chlorophyll a significantly decreased, water clarity significantly increased, and seagrass coverage in western Shinnecock Bay has expanded. Collectively, this study provides new insight into the manner in which bivalve spawner sanctuaries can be used to restore bivalve populations, improve water quality, and regrow seagrasses.

Evidence that Spawner Sanctuaries Caused Increased Landings and Densities of Hard Clams

Multiple lines of evidence demonstrate that the significant increase in the hard clam population in Shinnecock Bay can be attributed to the establishment of hard clam spawner sanctuaries. The more than 3.2 million hard clams planted in western Shinnecock Bay displayed undetectable mortality and patterns of gonadal condition indicative of annual spawning during the past decade (Doall et al., 2008). Physical circulation models demonstrated that larvae in the western part of the bay would be transported to eastern Shinnecock Bay within the two-week duration of the larval stage of M. mercenaria (Castagna, 2001) and also indicated that larvae spawned in eastern Shinnecock Bay would be tidally advected out to the Atlantic Ocean in less than a week. The quantification of hard clam larvae in Shinnecock Bay demonstrated that peak larval densities were found within the Tiana Bay spawner sanctuary at the time hard clams were spawning, and then appeared eastward by the Ponquogue Bridge three days later, and then in eastern Shinnecock Bay days after that, a pattern fully consistent with that shown by hydrodynamic models.

The timing of increases in recruitment and landings following the establishment of spawner sanctuaries also suggests that the spawner sanctuaries contributed to these increases. The planting of adult clams began in 2012 and rapidly continued through 2017, totaling 2 million clams planted by 2016 and 3 million by 2017. Increases in clam densities, particularly in the eastern Bay, were largely due to increases in small clams (<25 mm; < 3 years) produced during the time frame in which the spawner sanctuaries were established. Increases in commercial landings began to accelerate in 2016, four years after the initial establishment of the spawner sanctuaries and was driven primarily by an increase in landings of small clams. As previously noted, it takes ~three years for hard clams to reach harvestable size in waters on the south shore of Long Island (Kraeuter et al., 2005). While increased landings can simply be a function of increased harvest pressure, the differential changes in the size of clams landed (i.e., 26-fold increase in landings of littleneck clams versus 7-fold increase in landings of cherrystone and chowder sized clams) was also supportive of the hypothesis that the increase in harvested clams since 2012 originated from population increases caused by higher recruitment. And while the increased landings could be related to broad-scale environmental conditions favoring the expansion of hard clam populations in the region, neighboring lagoonal estuaries across the south shore of Long Island displayed unchanged or decreased landings during the same period landings increased 16-fold in Shinnecock Bay.

Differential Trajectories of Eastern and Western Shinnecock Bay

Despite the location of hard clam spawner sanctuaries in the western half of Shinnecock Bay and despite the prediction that hard clam larvae from these locations would circulate to both eastern and western Shinnecock Bay, the trajectories of hard clam populations in these two regions were different. While both halves of the Bay saw increases in the densities of hard clams since restoration began, there was a 18-fold increase in the eastern bay but only a 50% increase in the western bay where a portion of larvae were predicted to have settled. Survival rates of M. mercenaria larvae are very low (<1%; Fegley, 2001) as this and other early life stages of clams are highly sensitive to environmental conditions including acidification (Talmage and Gobler, 2009; Talmage and Gobler, 2010; Talmage and Gobler, 2011) and harmful brown tide blooms (Greenfield and Lonsdale, 2002; Gobler and Sunda, 2012). From 2007-2017, maximal and average cell densities of the harmful brown tide alga (1.7 x 10⁶ and 6.97 x 10⁴ cells mL^-1, respectively) in western Shinnecock Bay were 11- and 24-fold times higher, respectively, than in eastern Shinnecock Bay (1.48x 10⁵ and 2.88 x 10³ cells mL^-1; SCDHS, 1985-2021) likely due to the strong tidal flushing with Atlantic Ocean water in the eastern Bay. The absolute densities of brown tide are impactful as early life stage clam feeding slows at 3.5 x 10⁴ cells mL^-1 and ceases at 2 x 10⁵ cells mL^-1 (Bricelj et al., 2001), levels found in the west during the summers of 2012 – 2016 but never seen in the east (SCDHS, 1985-2021). Moreover, juvenile clams experience near complete mortality when exposed to > 10⁵ cells mL^-1 of brown tide algae for several weeks (Greenfield and Lonsdale, 2002), a level common in the west during summers from 2007-2016, but never seen in the east. Similarly, cruise tracks measuring pH during the clam spawning period (June-July) have mapped normal pH levels in the eastern Bay (7.9 – 8.1; 400 – 500 µatm) but low pH (< 7.5) levels in the western Bay (Figure S9), that can significantly increase the mortality of larval hard clams (Talmage and Gobler, 2009; Talmage and Gobler, 2010; Talmage and Gobler, 2011). In addition, recent field studies surveying the abundance of larval bivalves across various lagoonal sites revealed that the combination of high temperature and low DO reduced larval abundance (Weinstock and Collin, 2021), and western Shinnecock Bay is warmer and has lower DO than the eastern bay during the larval spawning season (SCDHS, 1985-2021). Further, benthic surveys suggest that densities of crabs are higher in the western Bay compared to the east (Figure S10). Hence, higher predation rates on juvenile clams (Juanes, 1992) as well as environmentally driven lower survival rates likely have both contributed to a slower recovery of the hard clam population in the western Bay compared to eastern Shinnecock Bay. It is possible that as hard clam population densities slowly increase within the western Bay and turnover times continue to decline, water quality will continue to improve, making environmental conditions more favorable for the survival of early life stage clams. It is also possible that organic-rich, acidic muds that minimize survival of newly set hard clam are an obstacle for recovery (Green et al., 2004; Green et al., 2009) in which case buffering of sediments with crushed shell may be needed to promote restoration of clam populations in some western regions (Curtin et al., 2022).

The Role of Biological Filtration in Facilitating Alternately Stable States

For decades, it has been recognized that bivalves are ecosystem engineers that, when present at high densities, can shape the character of estuaries by filter feeding (Jackson et al., 2001; Lotze et al., 2006) by transferring organic matter from the water column into sediments, increasing denitrification rates (Wozny, 2021), light penetration, and sediment nutrient stores to the benefit of seagrasses like Z. marina (Carroll et al., 2008; Wall et al., 2008) that in turn serve as a critical habitat for shellfish and juvenile fish (Heck and Orth, 1980; McDevitt-Irwin et al., 2016). During the end of the twentieth century, Shinnecock Bay was emblematic of an estuary that suffered the consequences of the loss of hard clam populations as landings had declined by 99.4% from 1970 to 2011 (NYSDEC, 1970-2011), densities were < 1 clams m^-2 (Weiss et al., 2007), dense brown tides occurred annually, and seagrass populations had been eradicated from much of the western Bay (NYSDEC, 2009). The strategic placement of hard clam spawner sanctuaries in regions shown to facilitate robust gametogenesis (Doall et al., 2008; Padilla and Doall, 2012) with suitable water and sediment quality and physical circulation patterns conducive for larval retention assisted in repopulating Shinnecock Bay with hard clams as well as significantly increasing the influence of their filter feeding on the estuary. The shallow nature of Shinnecock Bay (mean depth ~2 m) helps maximize the influence of bivalve filtration on the water column compared with deeper, more voluminous estuaries (Gerritsen et al., 1994; Pomeroy et al., 2006). While biological filtration times were three weeks and three months in western and eastern Shinnecock Bay, respectively, prior to the commencement of this project, the stocking and regrowth of hard clams across the estuary reduced the biological turnover time of the bay to 10 and 18 days, respectively, while the turnover times within spawner sanctuaries were less than three days. Prior studies have suggested a biological turnover time of three days or less was ideal for establishing biological control of algal blooms (Cerrato et al., 2004; Wall et al., 2008). While this was not achieved outside of the spawner sanctuaries, we note that tidal flushing likely achieves this standard in the eastern bay and acts additively with clam filtration rates and other mortality processes (e.g. other bivalves, zooplankton, etc.) in the western bay to provide removal rates that have contributed toward significant reductions in the levels of phytoplankton biomass in western Shinnecock Bay (p<0.0001), significant increases in water clarity (p<0.001), and the eradication of harmful brown tides from 2017-2021. This is the first time there has been five consecutive years without brown tide in any estuary on the south shore of Suffolk County, NY, since these events began in 1985 (SCDHS, 1985-2021). During the same time period in which water quality improved in western Shinnecock Bay, regions of the south shore estuary lagoon system in NY, USA, not influenced by spawner sanctuaries in Shinnecock Bay nor ocean inlets (i.e., central Great South Bay) saw densities of brown tide, levels of chlorophyll a, or Secchi disc depths increase in part to regional increases in nitrogen loads (Gobler et al., 2019). During the past decade there have been intensifying watershed nitrogen loads to south shore estuaries (SCSWWP, 2020) accounting for the declining water quality in Great South Bay and leaving increases in clam population densities in Shinnecock Bay as the most parsimonious explanation for improved water quality there.

These improvements in water clarity in western Shinnecock Bay were coincident with both the restoration of hard clam populations across the estuary and the regrowth of more than 4.74 x 10⁵ m² of Z. marina meadows from 2013-2017. This time frame is coincident with the time period when the clam-induced turnover time of the western bay declined from 21 to 11 days and suggests the ecosystem-wide impacts of this restored, dense population of filter feeding hard clams contributed to this change. Bio-optical modeling also revealed distinct increases in light penetration in western Shinnecock Bay, specifically within the regions of hard clam spawner sanctuaries where biological turnover times decreased to three days. Given the bio-optical model also demonstrated the water column light requirement for Z. marina has been extended to 2.4 x 10⁵ m² of bottom in western Shinnecock from 2017 to 2021 (data not shown), it is anticipated that the availability of clear aerial images of western Shinnecock Bay from 2021 may reveal a further expansion of seagrass during the period of the bio-optical model, 2017-2021.

Prospects for Estuarine Ecosystem Restoration

The restoration approach used for M. mercenaria in Shinnecock Bay was based on the premise that the population was recruitment limited and that improving fertilization success by densely co-planting millions of clams in individual sanctuaries with >25 individuals m^-2 would significantly improve reproductive output (Levitan, 1991; Mann and Evans, 2004). Similar approaches with other bivalves such as bay scallops (Argopecten irradians) have met with success in other regions of the US east coast (Peterson et al., 1996; Goldberg et al., 2000; Tettelbach et al., 2013; Tettelbach et al., 2015). The significant recovery of the hard clam population in Shinnecock Bay catalyzed the establishment of the Long Island Shellfish Restoration Program by the NYSDEC in 2018 using a highly similar approach to what has been reported here: creating small regions (~2,000m²) with high densities of hard clams (~25 m^-2) in no-harvest regions with water and sediment quality that is supportive of multiple life stages of hard clams and residence times that will allow larvae to settle in neighboring estuarine regions. Given the high mortality of juvenile hard clams (Juanes, 1992; Peterson et al., 1995; Gosselin and Qian, 1997; Munroe et al., 2015), it would seem that mass spawning events associated with spawner sanctuaries are required to overcome such mortality and is a preferrable restoration approach given the ability of each pair of adult hard clams to create millions of offspring annually (Kraeuter and Castagna, 2001). Moreover, the larger aerial coverage of larval settlement compared to that obtainable through manual planting of juvenile clams provides a greater likelihood of juvenile clams settling in regions hospitable for their survival. And while dense assemblages of bivalves in experimental and aquaculture settings that can lead to food limitation (Rheault and Rice, 1996; Wall et al., 2011), the planting densities used here were below such thresholds, a conclusion supported by the high survival and gametogenesis rates of clams within spawner sanctuaries.

Several attributes of the hard clam (M. mercenaria) made it an ideal choice for restoration of Shinnecock Bay. First, as a native species, there was a background, albeit diminished, population of individuals to build upon and there were not concerns associated with invasive species (e.g. zebra and quagga mussels in freshwater systems; Strayer, 2010). In addition, as a longer-lived bivalve (40+ yr), increases in M. mercenaria densities should be more stable than shorter-lived bivalves. For example, working with the shorter-lived bay scallop (Argopecten irradians; 2 yr), Tettelbach et al. (2013; Tettelbach et al., 2015) described a faster population recovery than was described here, but this shorter-lived population was less stable and has recently experienced sharp declines associated with climate change and other factors (Tomasetti, 2022).There are multiple attributes of shallow lagoons such as Shinnecock Bay that make them well-suited to the simultaneous restoration of bivalves and water quality. The shallow nature of this estuary (~ 2 m) maximizes benthic-pelagic coupling and, therefore, the ability of bivalves to exert maximal filtration pressure on this system as compared to deeper estuaries where such impact is minimized (Pomeroy et al., 2006). In addition, Shinnecock Bay is microtidal, with a typical tidal range of ~1 m, allowing bivalve filtration pressure to be maximal relative to tidal export of plankton (Banas et al., 2007) while also maximizing larval retention (Chen et al., 1997).Another key aspect of the repopulation of Shinnecock Bay with hard clams has been the use of spawner sanctuaries where clams cannot be harvested. No-take zones are a common and increasingly used restoration practice in marine ecosystems (Halpern, 2003; Lester and Halpern, 2008; Edgar et al., 2014; Gownaris et al., 2019; Grorud-Colvert et al., 2021). There have been multiple examples of hard clam populations being overfished in NY (Kraeuter et al., 2008) and elsewhere (Peterson, 2002). Given the extended time for clams to grow to harvestable size (at least three years; Kraeuter and Castagna, 2001), reductions in hard clam populations to recruitment limiting densities is a situation that may not be reversible without the introduction of new adult spawning stock (Peterson et al., 1996). The creation of spawner sanctuaries assures a continued supply of new recruits into ecosystems and, in some ways, mimics a situation present on Long Island in the middle of twentieth century when duck farms were prominent along south shore lagoonal estuaries (Great South Bay, Moriches Bay) and forced the closure of northern sections of bays to shellfishing due to fecal contamination (Ryther, 1954; Swanson et al., 2010). This allowed large adult clams to become established in closed regions that produced settlers that recruited and grew in the southern regions of these estuaries and contributed to the largest harvests of hard clams ever recorded in NY (McHugh, 1983). Through the latter half of the twentieth century, duck farms closed, fecal contamination was reduced, northern regions of estuaries were re-opened to shellfishing, and were subsequently overharvested (Swanson et al., 2010). With the loss of this spawning stock, the hard clam fisheries in these regions collapsed (Kraeuter et al., 2008). In some respects, the restoration approach used in Shinnecock Bay recreated the situation that supported the record setting harvests of the mid-twentieth centuries, with the difference being that today, spawner sanctuaries in Shinnecock Bay have been designated no-harvest zones for the purposes of restoration. This underscores the importance of political and stakeholder engagement, which is needed to establish and enforce such zones.

Longer Term Prospects and Perspectives

The longer life span of M. mercenaria (up to four decades; Kraeuter and Castagna, 2001) compared to other temperate bivalves (e.g., A. irradians, two years) makes it an ideal target for long-term restoration (Peterson, 2002). While this extended life span suggests M. mercenaria spawner sanctuaries could have a lasting impact on the Shinnecock Bay and similar coastal ecosystems, the rapid pace of climate change must be considered for restoration projects. The coastal waters of NY have been warming at two-to-three times global averages during summer since 1982 (Gobler et al., 2017; Griffith et al., 2019), and larval stages of M. mercenaria can experience enhanced mortality under elevated temperatures (28°C; Talmage and Gobler, 2011). This bivalve, however, has a geographic range that extends south into the Gulf of Mexico where individuals experience slowed growth associated with the warmest temperatures of the year (Kraeuter and Castagna, 2001). Hence, warming temperatures may alter the phenology of restored hard clam populations, with an increased prevalence of heat waves (IPCC, 2019) more likely to cause mortality in early life stage individuals (Talmage and Gobler, 2011; Stevens and Gobler, 2018) rather than simply progressive warming. Ocean acidification may be an even greater risk for restored hard clam populations as recent modeling estimated that the five-year risk of a 50% decline in M. mercenaria populations on the south shore of Long Island increases from 25% to 79% with an increase in pCO₂ levels from 400 μatm to 800 μatm (Grear et al., 2020). Revision of state and federally recommended water quality criteria to account for the multiple stressors associated with eutrophication and climate change may help improve the long-term viability of this fishery (Tomasetti and Gobler, 2020). Given that acidification in estuaries can be caused by eutrophication (Melzner et al., 2013; Wallace et al., 2014; Wallace and Gobler, 2021) regional efforts to combat excessive nutrient loading may also likely benefit hard clam spawner sanctuaries by lessening the intensity of harmful algal blooms (Heisler et al., 2008; Gobler et al., 2011) and the intensity of coastal acidification (Melzner et al., 2013; Wallace et al., 2014; Wallace and Gobler, 2021).

Conclusions

This study demonstrates the cascading, ecosystem-wide effects hard clam spawner sanctuaries can have on shallow lagoonal estuarine ecosystems. A decade of benthic and pelagic monitoring revealed that the stocking of no-harvest zones with dense assemblages of adult hard clams (27 m^-2) in areas with water and sediment quality supportive of gametogenesis and extended larval residence times caused a 18-fold increase in hard clam densities and a 16-fold increase in hard clam landings in Shinnecock Bay. Multiple lines of evidence linked the regrowth of clam populations to the establishment of spawner sanctuaries including the survival and spawning patterns of stocked clams, physical circulation patterns in the estuary, spatial-temporal patterns of spawning and clam larvae abundances, the timing of increases in recruitment and landings, the disproportionate increase in landings of smaller/younger clams, and differential landings within Shinnecock Bay compared to other regional estuaries. Higher clam densities caused biological filtration times of the bay to decrease from three weeks to 10 days in the western Bay contributing to concentrations of HABs and chlorophyll a significantly decreasing and water clarity and the extent of seagrass beds increasing. Given the trajectories of many bivalve populations across the globe (Jackson et al., 2001; Lotze et al., 2006; Beck et al., 2011; Ermgassen et al., 2012) and the recognition of the need for marine ecosystem restoration (Duarte et al., 2020), the approach taken here may be promising for recovering hard clam and other bivalve populations in estuaries around the globe. Finally, the use of protected sanctuaries in Shinnecock Bay and a long-lived bivalve species (M. mercenaria; 40+ yr) provides hope that effects of restoration presented here will persist over an extended period.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author Contributions

CG, MD, and BP conceived and designed the study. All authors contributed to the components of the study, including data collection, post-study analyses, analysis, and interpretation of the data. CG and BP contributed reagents, materials, and analysis tools. CG and MD wrote the paper. BP, AG, JC, CS, EP, and FD provided important edits. All authors contributed to the article and approved the submitted version.

Funding

This study was primarily support by the Laurie Landeau Foundation with matching funds from the Simons Foundation and additional support from the New York State Department of Environmental Conservation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We greatly appreciate the logistic support provided by the Stony Brook Southampton Marine Science Center staff throughout this study. We thank Jennifer O’Dwyer and Wade Carden from the New York State Department of Environmental Conservation for providing the historical hard clam landings data for New York State.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.911731/full#supplementary-material

References