Agricultural nitrogen (N) contributes a dominant percentage to global N pollution in the coastal zone. Emerging research on N isotopes in bivalve shells has shown value for reconstructing historical increases in estuarine wastewater inputs. However, applications for fertilizer N are understudied. Here, we integrate the study of organic N isotopes, in concert with δ¹⁸O and δ¹³C, in estuarine bivalve shells to investigate spatial and long-term change in nitrogen inputs and sources. Modern, museum collection, and subfossil specimens of the genera Mytilus and Ostrea were profiled in a California estuary with an intensely agricultural watershed. Spatial patterns in bivalve isotopic composition reflected gradients in watershed nutrient inputs and productivity parameters. Furthermore, comparison of modern and historical periods revealed changes in nutrient source or processing over the last 1000 years. The N isotope values from shells offer perspective on agricultural pollution in estuaries.

Elkhorn Slough (36° 48' 25.18"N, 121° 47' 23.24"W) is a shallow estuary (2.5 m deep) that extends inland 10 km from Monterey Bay in central California, USA. Elkhorn Slough drains a watershed of 585 km², and while a portion has been protected, 24% of the watershed is cultivated, generally with crops that have moderate to high fertilizer demands and/or are multi-cropped. Outside of tidally restricted areas or during infrequent winter storms, the salinity and temperature of the Elkhorn Slough’s estuarine waters are similar to those found in Monterey Bay (S=30 PSU; T=14^oC). It serves as a model system for studying agricultural N inputs. Water sampling has frequently documented high nitrate levels exceeding 1000  μM, an order of magnitude greater than many other well-studied eutrophic estuaries. Furthermore, historical changes in the hydrology of the estuary, including creation of an artificial oceanic inlet at the harbor in 1946, may have also altered marine nutrient exchange.

In October of 2021, live Mytilus specimens were collected at 11 sites throughout the estuary, while recently deceased Ostrea lurida shells were collected from one location near the mouth of the estuary. Mytilus specimens were not identified to the species level because Mytilus galloprovincialis and Mytilus trossulus species hybridize in coastal California. Despite the scarcity of Ostrea, they were sampled due to extensive prior research on δ¹⁵N_shell values in oysters. Five specimens were collected at each site, based on prior research suggesting this as an optimal sample size for reproducibility.

Modern bivalve isotopes were related to water chemistry measured in grab samples collected monthly from November 2020 to October 2021 by the National Estuarine Research Reserve. Several parameters were measured in situ, including temperature, dissolved oxygen, and salinity, using EXO2 multiparameter sondes (YSI, Yellow Springs, OH, USA). Grab samples were stored on ice, filtered within 48 hours, and analyzed for chlorophyll a using fluorometry and dissolved inorganic nitrogen using a Lachat QC8000 Flow Injection Analyzer at Moss Landing Marine Laboratories, following standard methods.

Historical specimens of Mytilus and Ostrea species (collected between 1933 and 1999) were obtained from the California Academy of Sciences museum. Fragments of Mytilus shells (dating from approximately 964 to 1760) were extracted from three archived sediment cores collected in 2010, which were previously measured for sedimentary isotopes on the top 50 cm. Chronological control for subfossil specimens was established using AMS radiocarbon dating of organic material (seeds and peat), generating a composite chronology from seven cores collected throughout the estuary to mitigate complications from the marine reservoir effect. The chronology was developed using a Bayesian age model called “Plum” that combines lead-210, radiocesium, and radiocarbon dating.

We selected moderately sized Mytilus specimens (~40mm) to control for age impacts on increasing δ¹⁵N_shell values. We sampled a narrow portion of recent calcite growth (average 4.2mm) at the ventral margin of the valve from the outer shell layer (OSL) of each specimen. Based on estimated growth rates, this sample size corresponds to approximately one year of growth. Ostrea shells were sampled from the inner shell surface. For both species, sampling was conducted on calcite to control for fractionation between mineral forms. Nitrogen content of the historical specimens was compared to modern samples to investigate degradation of organic matter over time.

Shells were sampled using a dental drill with a 1 mm tungsten carbide conical bit at a low speed of 1000 rotations min^-1 to minimize isotopic fractionation. Subfossil Mytilus material preserved in sediment cores was sampled using the same approach as modern samples, grinding 4-5 mm portions along the margin where preserved, or the full sample for fragments less than 5 mm long.

Samples (300-400 mg) were analyzed for δ¹⁸O_shell and inorganic δ¹³C_shell using an Elementar Isotope Ratio Mass Spectrometer (IRMS) (Elementar Americas Isoprime100) coupled with a Multiflow carbonate prep system. Carbonates were digested with 103% orthophosphoric acid at 70°C for four hours. The same samples (20-50 mg) were analyzed for δ¹⁵N_shell on the Elementar IRMS coupled with a Pyrocube elemental analyzer equipped with a carbon trap (NaOH) between the reduction and absorption tubes to remove excessive CO₂ and reduce interference of the carbonate matrix with organic isotope analysis. Standard materials represented 20% of the samples (run precision ± 0.19 ‰ for δ¹³C_shell, ± 0.41 ‰ for δ¹⁸O_shell, ± 0.71 ‰ for δ¹⁵N_shell). For inorganic samples, internal standards were calibrated against the international standards NBS-18 and NBS-19, achieving a precision of ± 0.09‰ (C) and ± 0.16‰ (O) based on long-term replication. For organic samples, calibration was performed with IAEA N1, IAEA N3, and USGS 26 standards, with precision of ± 0.32 ‰. Internal standard “Elk” (1.05% N) closely matches the N content of the Mytilus and Ostrea samples. Every tenth sample was duplicated. We also measured δ¹⁵N _tissue (freeze-dried and ground soft tissues of Mytilus), δ¹⁵N_SPM (suspended particulate matter), and δ¹⁵N_NO3 (water column nitrate) values at several select sites to verify the δ¹⁵N_shell patterns (Fig. S8). Isotopic data are expressed in parts per thousand (‰) using the delta notation (δ¹⁵N, δ¹⁸O, and δ¹³C), where δ = [(R_sample−R_standard)/R_standard] × 1000.

We used a Partial Least Squares (PLS) model in R version 4.3.1 to relate modern δ¹⁵N_shell values with spatial patterns of modern water chemistry (‘plsr’ package). This modeling approach was chosen for its ability to handle multicollinearity, which is typically present in water chemistry data. The isotopic signatures were compared to one year of water chemistry data, represented as median values, which are temporally comparable to amount of shell material sampled. Predictors of δ¹⁵N_shell values included δ¹⁸O_shell, δ¹³C_shell, NO₃, NH₄⁺, temperature, salinity, dissolved oxygen, turbidity, and distance to the oceanic inlet. Statistical analysis of change over time was performed by grouping data into bins according to time: (1) modern samples collected in 2021; (2) museum specimens from 1933 − 1999 C.E.; and (3) core specimens from circa 964 to 1760 C.E. The time bins were compared using Analysis of Variance (ANOVA) because the data are normally distributed, as confirmed by the Shapiro–Wilk Test and histograms.