Seasonal to Inter-Annual Variability of Primary Production in Chesapeake Bay: Prospects to Reverse Eutrophication and Change Trophic Classification

Estuarine-coastal ecosystems are rich areas of the global ocean with elevated rates of organic matter production supporting major fisheries. Net and gross primary production (NPP, GPP) are essential properties of these ecosystems, characterized by high spatial, seasonal, and inter-annual variability associated with climatic effects on hydrology. Over 20 years ago, Nixon defined the trophic classification of marine ecosystems based on annual phytoplankton primary production (APPP), with categories ranging from “oligotrophic” to “hypertrophic”. Source data consisting of shipboard measurements of NPP and GPP from 1982 to 2004 for Chesapeake Bay in the mid-Atlantic region of the United States supported estimates of APPP from 300 to 500 g C m−2 yr−1, corresponding to “eutrophic” to “hypertrophic” categories. Here, we developed generalized additive models (GAM) to interpolate the limited spatio-temporal resolution of source data. Principal goals were: (1) to develop predictive models of NPP and GPP calibrated to source data (1982 to 2004); (2) to apply the models to historical (1960s, 1970s) and monitoring (1985 to 2015) data with adjustments for nutrient loadings and climatic effects; (3) to estimate APPP from model predictions of NPP; (4) to test effects of simulated reductions of phytoplankton biomass or nutrient loadings on trophic classification based on APPP. Simulated 40% decreases of euphotic-layer chl-a or TN and NO2 + NO3 loadings led to decreasing APPP sufficient to change trophic classification from “eutrophic’ to “mesotrophic” for oligohaline (OH) and polyhaline (PH) salinity zones, and from “hypertrophic” to “eutrophic” for the mesohaline (MH) salinity zone of the bay. These findings show that improved water quality is attainable with sustained reversal of nutrient over-enrichment sufficient to decrease phytoplankton biomass and APPP.

www.nature.com/scientificreports www.nature.com/scientificreports/ of phytoplankton biomass or nutrient loadings on trophic classification based on APPP. Targeting these goals allowed us to test the hypothesis that nutrient reductions adopted by Chesapeake Bay jurisdictions would lead to reduced APPP and a change of trophic classification.

Methods
cruises. Data were collected on 78 cruises on four research vessels from March 1982 to November 2004. All stations we occupied to measure NPP and GPP are depicted in Fig. 1. Data from a subset of these cruises were analyzed previously 7 . Experimental protocols for measuring NPP and GPP were specific to individual projects identified by acronyms in Table 1. CB cruises  consisted of an initial transect for horizontal mapping to determine salinity, temperature, chl-a, nutrient, and turbidity gradients and to select stations for measuring NPP. ProPhot and FITS cruises  and LMER PROTEUS cruises (1989-94) occupied stations along a north-south transect to map water-quality properties and to measure NPP and ancillary properties. NASA cruises (1993-94) occupied three stations per cruise in the lower bay, plume, and adjacent shelf waters interspersed with mapping legs to measure NPP, GPP, and ancillary properties. LMER TIES cruises (1995)(1996)(1997)(1998)(1999)(2000) sampled stations at 0.5° latitudinal increments, and additional stations lateral to the north-south axis, for a total of nine to 27 stations per cruise. LMER TIES cruises also included legs for horizontal and vertical mapping, in-situ sampling of  www.nature.com/scientificreports www.nature.com/scientificreports/ Vertical profiles. Vertical profiles of salinity, temperature, dissolved oxygen (DO), and chl-a fluorescence were determined from hydro-casts at a set of stations spaced 10-20 km along horizontal transects. Some of the methods to obtain vertical profiles of these properties were described previously 7 . Hydrocasts were made with Yellow Springs salinity and DO meters (see above) on 1982-1983 cruises, a Sea Bird model 9 CTDFO 2 on 1987-1988 cruises, and a Neil Brown Mark III CTDFO 2 on 1989-2004 cruises. Discrete samples were collected in Niskin bottles to calibrate fluorometers and DO meters. Surface sunlight during the day (E 0 = downwelling irradiance, photosynthetically available radiation = PAR) was measured with a Li-Cor quantum meter model 190 S (or equivalent) coupled to a Li-Cor model 550 or 1000 integrator. The sensor was mounted near the deck incubators we used to measure simulated in-situ NPP and GPP. Diffuse light attenuation coefficient for PAR (K PAR ) was determined from vertical profiles with a submersible Li-Cor quantum meter model 188B with a 192S sensor or equivalent for CB, ProPhot, FITS, and LMER PROTEUS, TIES, and NASA cruises. Secchi depths (m) were determined for all stations and cruises. Euphotic-layer depth (Z p ) was estimated as the depth to which 1% of E 0 penetrated based on vertical profiles of downwelling irradiance (E z ), or from Secchi depth using calibration regressions 32 . NSF LMER TIES and NASA cruises also conducted vertical profiles with a Biospherical Instruments multi-spectral environmental radiometer (MER-2040/41) to measure K PAR and spectral light attenuation.
Water-quality properties. Chl-a was determined using spectrophotometric and fluorometric measurements on acetone extracts (90%) of particulate material collected by vacuum filtration onto glass-fiber filters (Whatman GF/F or equivalent, 0.3-0.8 µm nominal pore sizes). Spectrophotometric chl-a was derived from trichromatic equations applied to absorbances measured on a Beckman DK-2 or equivalent, and fluorometric chl-a was measured on a Turner model 110, 111, or Turner Designs model 10 calibrated by spectrophotometry 20,21 . Secchi depth was the depth at which a 30-cm white disk became invisible when lowered over the side of the research vessel. NO 2 + NO 3 was measured using analytical methods documented by the EPA/CBP 33,34 following protocols described by D'Elia et al. 35 .
Simulated in-situ incubations. NPP was measured at 723 stations from 1982 to 2004, and GPP at 525 stations from 1995 to 2004. Whole-water samples were collected in Niskin bottles mounted on a rosette sampler at sunrise at a depth of 0.5 to 1.0 m, contents were pooled in a darkened carboy, and dispensed to 125-150 ml glass incubation bottles. The euphotic layer is well mixed in the bay and chl-a is homogenously distributed in the upper 5-10 m. NPP and GPP were determined using 14 C-sodium bicarbonate uptake in deck incubators cooled with flowing surface water (±1 °C of in-situ temperatures). 2 to 5 µCuries of 14 C-sodium bicarbonate (ICN Pharmaceuticals, Inc., or Amersham Searle, Inc.) were added to each incubation bottle. Total 14 C-sodium bicarbonate activity was determined for a time-zero aliquot from one of the incubation bottles, and from a small amount of stock isotope added to scintillation cocktail made basic with NaOH (Aquasol, New England Nuclear, Inc., or equivalent). Incubation bottles were exposed to a range of sunlight levels using neutral density screens providing 100 (no screens), 58, 34, 21, 11, 4 and 1% transmission to simulate seven depths in the euphotic layer. Dark uptake was measured in an opaque bottle and used as a tare value to adjust uptake in illuminated bottles.
NPP on cruises from 1982 to 1983 was measured in 3-4 h incubations repeated four times during the photoperiod, with additional bottles incubated throughout the night. NPP on cruises from 1987 to 2004 was measured in 24-h incubations, and GPP on cruises from 1995 to 2004 was measured in 4-6 h incubations. Concurrent measurements using 14 C and O 2 methodologies confirmed this approach accurately estimated NPP and GPP 7 . Duplicate 25-150 ml subsamples (depending on phytoplankton biomass) were withdrawn at the end of these periods and filtered onto glass fiber filters (Gelman AE or Whatman GF/F) under low vacuum pressure (<150 mm Hg). Filter pads were rinsed with filtered water of equivalent salinity as the sample and acidified with 0.01 N HCl in a fume hood to remove residual inorganic label. Activities were determined on a Packard Instruments Tri-Carb or model 3320 liquid scintillation counter. Duplicate aliquots were withdrawn from incubation bottles to determine chl-a using methods described above. Total CO 2 was measured by gas-stripping, capture, and analysis on a Beckman model 864 infrared analyzer for 1982-83 cruises, Gran titration for 1987-1988 cruises, and by gas chromatography on a Hach Cable Series 100 AGC for 1989-2004 cruises.
computations of npp and Gpp. Equations to compute NPP and GPP combined terms for 14 C-uptake in lighted bottles tared against non-photosynthetic 14 C-uptake in dark bottles, divided by the total 14 C activity added, and multiplied by a discrimination factor of 1.05 for 14 C vs 12 C and a term for total CO 2 (mM) converted to weight (mg m −3 ). The resulting volumetric rates (PP, mg C m −3 h −1 ) were used to compute NPP and GPP by converting simulated incubation depths as percent surface irradiance (E 0 ) to actual depths based on K PAR from vertical profiles of irradiance (E d ) or Secchi depth. Multiple-segment trapezoidal integration was applied to PP from the surface to euphotic-layer depth to obtain NPP and GPP. NPP was determined from 24-h incubations and GPP from 4-6 h incubations scaled to the photoperiod using E 0 during incubations as a proportion of total E 0 for the day. Optimal photosynthesis, P B opt , was determined as maximum PP in simulated in-situ incubations normalized to chl-a. Observed values of log 10 P B opt were binned in 1° increments as a function of SST following the approach of Son et al. 31 and analyzed by polynomial regression. Estimates of log 10 P B opt from these regressions were used as model inputs to predict NPP for historical (1960s, 1970s) and monitoring  programs.
Freshwater discharge, climatic effects. Daily freshwater discharge from the Susquehanna River (SRF) and monthly cumulative discharge (SUM) were obtained from the U.S. Geological Survey (USGS) (http://md.water.usgs.gov) for the Conowingo Dam gaging station near the head of the bay (latitude 39° 39′ 28.4″, longitude 76° 10′ 28.0″). TN and NO 2 + NO 3 loadings were obtained as monthly values (10 6 kg) from the USGS. We focused on the Susquehanna River as this large river dominates distributions of nutrients and phytoplankton in the main stem bay. Nutrient loadings from other tributaries are significantly reduced by processes within their confines, (2020) 10:2019 | https://doi.org/10.1038/s41598-020-58702-3 www.nature.com/scientificreports www.nature.com/scientificreports/ attenuating effects of lateral inputs on the bay proper [16][17][18] . Climatic effects were quantified by applying GAM to input files for NPP, GPP, and water-quality properties, including terms for freshwater flow (log 10 monthly SRF, log 10 monthly SUM) and salinity. Model predictions of NPP and GPP for low-flow, "dry" conditions were based on flow terms set at 10 th percentiles joined by salinity terms at 90 th percentiles; mean-flow model predictions were based on flow and salinity terms held constant at their mean values; model predictions for high-flow, "wet" conditions were based on flow terms set at 90 th percentiles joined by salinity terms at 10 th percentiles.

Analytical steps.
A flow-chart of analytical steps, including data sources, modeling approach, and simulations is presented in Fig. 3. Statistical analyses were conducted using "R" v. 3.6.0. Non-linear fits of time-series data were developed using GAM from the "R" package 'mgcv' , containing functions similar to those designed by T. Hastie in S-Plus, and based on a penalized regression-spline approach including automatic smoothness selection [36][37][38] . Model predictions of NPP and GPP were adjusted for climatic effects using GAM, as described previously for water-quality properties [16][17][18][19] . A comparable approach by Beck and Murphy 39 compared GAM to weighted regressions of time, discharge and season (WRTDS) developed by Hirsch et al. 40 , noting the flexibility of GAM to add relevant predictor variables as a strength of GAM over alternatives.
Model fits, residuals, flow-adjusted model predictions at monthly increments, adjusted R 2 , generalized cross validation (GCV) score, % deviance explained, p-values for F-statistics, and root mean square error (RMSE) were obtained for each model. Degrees of smoothing (knots = k) were selected by the "R" package 'mgcv' to minimize the GCV score, followed by post-hoc adjustments of "k" for individual terms using the function "gam. check". Graphical presentations were prepared with Kaleidagraph 4.5.2 (Synergy Software, Inc.). These include time series of mean, monthly NPP and euphotic-layer chl-a ( Fig. 4a-c), polynomial regressions of P B opt on SST ( Fig. 5a,b), observed vs model-fitted values of log 10 NPP and log 10 GPP ( Fig. 6a-f), probability distributions of observed and predicted log 10 NPP and log 10 GPP ( Fig. 7a-d), comparisons of predicted log 10 NPP from gam1 and gam2 (Fig. 8a-c), time-series of log 10 chl-a, log 10 euphotic-layer chl-a, and model predictions of log 10 NPP from 1985 to 2015 ( Fig. 9a-i), flow-adjusted model predictions of log 10 euphotic-layer chl-a and log 10 NPP for low-flow, mean-flow, and high-flow conditions ( Fig. 10a-f), time series of APPP based on mean-flow model predictions of NPP and simulated reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings, and observed euphotic-layer chl-a ( Fig. 11a-c), and historical reconstructions of mean, monthly log 10 NPP and log 10 chl-a for the 1960s and 1970s with estimates of APPP ( Fig. 12a-f).

Results
Annual cycles. Spatial and seasonal differences were evident in mean, monthly euphotic-layer chl-a and NPP based on shipboard measurements in Chesapeake Bay from 1982 to 2004 (Fig. 4a-c). Seasonal maxima of phytoplankton biomass and production occurred in summer for the OH salinity zone, with euphotic-layer chl-a ~40 mg m −2 and NPP ~1200 mg C m −2 d −1 in July and August (Fig. 4a). Contrasting patterns for the MH salinity zone consisted of a well-developed spring bloom with euphotic-layer chl-a >80-100 mg m −2 in April and May, displaced several months from a summer maximum of NPP >2000 mg m −2 d −1 in July (Fig. 4b). A second maximum of euphotic-layer chl-a from 60-80 mg m −2 occurred in fall for the MH salinity zone but was not matched by a maximum of NPP (Fig. 4b). Annual cycles of euphotic-layer chl-a and NPP for the PH salinity zone were similar in profile and somewhat lower compared to those for the MH salinity zone, with maxima of euphotic-layer chl-a ~60 mg m −2 and NPP ~1700 mg m −2 d −1 in May and September, respectively (Fig. 4c). Table 2 summarizes the statistical properties of shipboard measurements from 1982 to 2004 by salinity zone and season, including Z p , salinity, SST, surface chl-a, euphotic-layer chl-a, NPP, and GPP. www.nature.com/scientificreports www.nature.com/scientificreports/ Models of p B opt , npp, and Gpp. P B opt was estimated from third-order polynomial regressions of observed log 10 P B opt on binned SST as described by Son et al. 31 . These regressions had p < 0.001 and R 2 > 0.80 (Fig. 5a,b). Estimates of log 10 P B opt from these regressions were combined with input data on water-quality properties to predict NPP and GPP. Table 3 presents a complete list of predictor variables for models of NPP and GPP. Simple, linear regressions of observed vs model-fitted log 10 NPP and GPP had p < 0.001 for OH, MH, and PH salinity zones (Table 4; Fig. 6a-f). Probability distributions of observed and model-fitted log 10 NPP and GPP confirmed that the models generated unbiased estimates as model predictions displayed statistical attributes indistinguishable from observations ( Fig. 7a-d).
Model predictions. Five GAM formulations were developed to accommodate estimates of log 10 NPP and log 10 GPP based on availability of data for predictor variables in different time periods (Table 5). We focused primarily on model predictions of log 10 NPP as these supported estimates of APPP. Model predictions of log 10 NPP were compared for gam1 and gam2, i.e., model formulations with and without the predictor variable "sequential year" (Seq_year). By omitting Seq_year in gam2, we avoided an assumption that long-term trends in calibration data (1982 to 2004) were unchanged for water-quality properties outside that range (2005 to 2015) that were used to predict log 10 NPP. Both gam1 and gam2 captured seasonal to interannual variability of log 10 NPP, with differences in long-term trends expressed in simple, linear regressions of model predictions on year (Fig. 8a-c). Subsequent analyses predicted log 10 NPP from 1985 to 2015 using gam2 based on water-quality properties as data inputs. An alternate model to gam2 was developed as gam3 to address data limitations for historical periods, using the term log 10 chl-a for biomass in place of log 10 euphotic-layer chl-a. gam3 was further modified to gam4 by omitting the term K PAR and to gam5 by omitting the term for season. These models were applied to input data from the 1960s and 1970s as log 10 euphotic-layer chl-a, K PAR , or season were either absent or too sparse to support model predictions with reasonable sample sizes. Climatic effects. Climatic effects on log 10 euphotic-layer chl-a and log 10 NPP are expressed as model predictions for low-flow, "dry", mean-flow, and high-flow, "wet" conditions for OH, MH, and PH salinity zones from 1985 to 2015 (Fig. 10a-f). Simple, linear regressions showed consistent upward trends of log 10 euphotic-layer chl-a for all three salinity zones (Fig. 10a-c). Low-flow conditions did not affect log 10 euphotic-layer chl-a for the OH salinity zone but led to decreased log 10 euphotic-layer chl-a for MH and PH salinity zones. High-flow conditions led to decreased log 10 euphotic-layer chl-a for the OH salinity zone, did not affect log 10 euphotic-layer chl-a for the MH salinity zone, and led to increased log 10 euphotic-layer chl-a for the PH salinity zone. NPP was less sensitive to climatic effects than euphotic-layer chl-a, documented as time series of log 10 NPP for low-flow, mean-flow, and high-flow conditions (Fig. 10d-f). Simple, linear regressions of mean-flow predictions of log 10 NPP showed upward trends for OH and MH salinity zones, but no trend for the PH salinity zone. euphotic-layer chl-a, Appp. Time series of mean, annual euphotic-layer chl-a based on observations from 1985 to 2015 showed upward trends for OH and MH salinity zones, but not for the PH salinity zone (Fig. 11a-c). Corresponding APPP estimates from mean-flow model predictions of NPP showed upward trends from 1985 to 2015 for all three salinity zones. APPP ranged from ~200 to 300 g C m −2 yr −1 for the OH salinity zone, ~250 to 550 g C m −2 yr −1 for the MH salinity zone, and 280 to 350 g C m −2 yr −1 for the PH salinity zone. Using Nixon's trophic classification 3 , APPP corresponded to "eutrophic" conditions for OH and PH salinity zones, and "hypertrophic" conditions for the MH salinity zone. APPP based on mean-flow model predictions of NPP with simulated 40% reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings showed reductions for all three salinity zones (Fig. 11a-c). These reductions of biomass or nutrient loadings changed APPP to mesotrophic conditions for OH and PH salinity zones, and to eutrophic conditions for the MH salinity zone. www.nature.com/scientificreports www.nature.com/scientificreports/ Historical reconstructions. Archival data on water-quality properties, including P B opt from polynomial regressions on SST (Fig. 5a), log 10 chl-a, and other predictor variables (Table 3), supported model predictions of log 10 NPP in the 1960s and 1970s (Fig. 12a-f). Model predictions from gam4 supported historical reconstructions of log 10 NPP for low-flow, "dry", mean-flow, and high-flow, "wet" conditions to capture climatic effects (Table 5). APPP based on mean-flow predictions of NPP for the OH salinity zone reached 555 g C m −2 y −1 in the 1960s, compared to 288 g C m −2 y −1 in the 1970s (Fig. 12a,d). APPP for MH and PH salinity zones showed similar patterns, ranging from 392 to 504 g C m −2 y −1 in the 1960s, and from 132 to 334 g C m −2 y −1 in the 1970s (Fig. 12b,c,e,f). Mean-flow predictions of log 10 NPP in the 1960s showed summer maxima for OH and MH salinity zones (Fig. 12a,b), while sparse data for the PH salinity zone limited resolution for that period (Fig. 12c). Observed log 10 chl-a showed summer maxima for OH, MH, and PH salinity zones, and an absence of a spring maximum for the MH salinity zone (Fig. 12b).

Discussion
npp and Gpp models. Important advances in NPP and GPP models for estuarine-coastal waters were made possible by increasingly sophisticated approaches and availability of calibration data. Common elements of production models dating to the 1950s include terms for photosynthetic efficiency, phytoplankton biomass, and light availability [41][42][43][44][45][46][47][48][49][50][51][52] . Light-utilization models specific to estuarine-coastal ecosystems, such as Chesapeake Bay, San Stepwise regressions of log-transformed variables from VGPM led to CBPM that supported estimates of NPP from aircraft and satellite ocean-color data 31,53 .
Here, we departed from VGPM and CBPM to develop production models for the bay using GAM. Selection of predictor variables for NPP and GPP models was informed by earlier studies on a subset of these data from LMER TIES cruises ( Table 1). Analysis of variance (ANOVA) showed that season and region explained most of the variability of phytoplankton properties, including NPP, chl-a, and floral composition 54 . Long-term (six-year) means of NPP were negatively correlated with the fraction of chl-a in diatoms, a property stimulated by high-flow, "wet" conditions. Multiple linear regression and principal component analysis identified SRF as a 'master variable' driving inter-annual variability of these properties. An important advantage of transitioning to GAM was the flexibility to include predictor variables such as SRF to adjust for climatic effects on hydrology and distinguish variability from trends.
Model forms developed here were guided by these findings, leading us to include predictor variables for salinity zone, salinity, month, season, year, and flow terms log 10 SRF and log 10 SUM (Table 3). These models proved effective to estimate NPP and GPP, exemplified by simple, linear regressions of observed vs. modeled log 10 NPP from gam2 with R 2 > 0.96 (Fig. 6a-c), and log 10 GPP from gam2 with R 2 > 0.80 (Fig. 6d-f). These fits are comparable to NPP estimates using CBPM with measured values of P B opt , and GPP estimates using CBPM with estimated values of P B opt 7 . Predictor variables in models of NPP for OH, MH, and PH salinity zones consistently showed highest F-values and lowest p-values for log 10 P B opt and log 10 euphotic-layer chl-a (Table 3). Several other terms were also significant predictor variables, i.e., TN and NO 2 + NO 3 loadings (OH salinity zone), salinity, month (MH salinity zone), and salinity, K PAR , SRF, and month (PH salinity zone).
Climatic effects. Our group has focused on climatic effects on hydrology impacting water quality and phytoplankton in recent studies of Chesapeake Bay [16][17][18][19] . Adolf et al. 54 explored this theme previously, reporting predictable consequences of SRF on phytoplankton dynamics. Statistical models based on long-term data extended these findings, documenting climatic effects on chl-a, floral composition, and NPP [16][17][18][19] . A logical sequence emerged from these studies wherein seasonal to interannual variability of freshwater flow and N loading regulates spatio-temporal distributions of phytoplankton [16][17][18][19] , consistent with the conclusion by Malone et al. 22 that P plays a limited, transient role in the OH salinity zone of the bay, while N limits phytoplankton biomass and production on the ecosystem scale. www.nature.com/scientificreports www.nature.com/scientificreports/ Despite evidence from shipboard, aircraft, and satellite data linking freshwater flow to phytoplankton dynamics in land-margin ecosystems, previous NPP and GPP models did not contain explicit terms for climatic effects on hydrology 12-14,51,52,54-60 . Analyses described here addressed this shortcoming, based on observations in Chesapeake Bay spanning several decades. Specifically, low-flow, "dry" conditions produce a landward shift of N-limitation toward OH and MH salinity zones, lower chl-a, lower NPP, and a lower proportion of diatoms in the phytoplankton flora; high-flow, "wet" conditions extend the area of N sufficiency seaward to MH and PH salinity zones, leading to higher chl-a, higher NPP, and a higher proportion of diatoms [16][17][18][19]54,60 . Climatic effects on bio-optically active constituents similarly affect light-limitation as higher inputs of dissolved and suspended materials occur for high-flow, "wet" conditions than for low-flow, "dry" conditions 16,18 . This latter observation may contribute to lower sensitivity of NPP than chl-a to climatic variability reported here (Fig. 10a-f).
Development of numerical water-quality criteria followed this logic, leading to model predictions that distinguished long-term trends from spatio-temporal variability 18 . Freshwater flow from the Susquehanna River, and frequencies of predominant weather patterns defined "dry" and "wet" conditions 53,[60][61][62] , and statistical models conditioned on specific input terms for flow and salinity supported predictions of mean, monthly chl-a, Secchi depth, and NO 2 + NO 3 19 . Here, we extended this approach to NPP and GPP models by including terms to adjust Figure 8. (a-c) Comparisons of log 10 NPP using predictions from generalized additive models (GAM) from 1985 to 2015; gam1 included an predictor variable for the time term "Seq-year" and gam2 omitted this variable (see Table 4). (2020) 10:2019 | https://doi.org/10.1038/s41598-020-58702-3 www.nature.com/scientificreports www.nature.com/scientificreports/ for climatic effects on hydrology (Tables 3, 5; Fig. 10a-f). This approach benefited from the flexibility of GAM to incorporate predictor variables traditionally used in production models, i.e., P B opt , chl-a or euphotic-layer chl-a, Z p , and SST, and to add variables for salinity zone, salinity, season, SRF, and TN and NO 2 + NO 3 loadings.
Appp. Cloern et al. 4 published a synthesis of APPP for estuarine-coastal ecosystems based on a comprehensive survey of the scientific literature. APPP for 131 ecosystems ranged from 105 to 1890 g C m −2 yr −1 , with a mean of 252 g C m −2 yr −1 . Ten-fold variability occurred within ecosystems and five-fold from year to year, with only eight time-series covering longer than a decade. One of the best-studied ecosystems in the survey was the Rhode River, a small sub-estuary adjacent to the MH salinity zone of Chesapeake Bay. Long-term measurements of photosynthesis by Gallegos 63 supported estimates of APPP ranging from 152 to 612 g C m −1 y −1 in the Rhode River, with a mean of 328 g C m −1 y −1 . APPP maxima occurred in years with dense spring blooms of Prorocentrum cordatum (formerly P. minimum) a dinoflagellate species that commonly forms "mahogany tides". Complex interactions of local and remote nutrient inputs affected the relationship of APPP in the Rhode River to SRF. High-flow conditions displaced the turbidity maximum, usually located in the OH salinity zone, south of the Rhode River mouth, causing elevated turbidity in the sub-estuary, washout of phytoplankton, suppression of the spring bloom, and decreased APPP.
Models of NPP calibrated with long-term measurements in Chesapeake Bay from 1982 to 2004 supported multi-year estimates of APPP, based on water-quality properties from 1985 to 2015 as model inputs. These estimates of APPP allowed us to resolve inter-annual variability for a three-decade span, rarely possible for Nixon's trophic classification, historical context. Mean-flow predictions of NPP were used to estimate APPP from 130 to >600 g C m −2 y −1 for OH, MH, and PH salinity zones, with increases from 1985 to 2015 matching euphotic-layer chl-a (Fig. 11a-c). APPP of this magnitude corresponds to "eutrophic" for OH and PH salinity zones, and "hypertrophic" for the MH salinity zone using Nixon's 3 trophic classification (Fig. 1). According to Cloern et al. 4 , Chesapeake Bay ranks among estuarine-coastal ecosystems that are heavily impacted by nutrient over-enrichment (their Fig. 4). We evaluated prospects for changing trophic classification based on APPP by simulating 40% reductions of biomass or nutrient loadings in models of NPP. These reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings led to decreased APPP and changed trophic status from "hypertrophic" www.nature.com/scientificreports www.nature.com/scientificreports/ to "eutrophic" for the MH salinity zone, and from "eutrophic" to "mesotrophic" for OH and PH salinity zones (Fig. 11a-c).
Simulated 40% reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings were based on goals established by the 1987 Chesapeake Bay Agreement 64 to reduce phytoplankton biomass sufficiently to reverse summer anoxia. Several interventions by management began in the 1980s when states bordering the bay banned phosphate in laundry detergents. Subsequent nutrient-management legislation was adopted by Maryland, Virginia, and Pennsylvania in the 1990s, aimed at reducing the over-application of commercial fertilizers and manure on agricultural lands. In 2004, the six states in the watershed, the District of Columbia, and U.S. EPA reached an agreement on comprehensive wastewater treatment permits, leading to numerical annual loading limits for over www.nature.com/scientificreports www.nature.com/scientificreports/ 470 municipal and industrial wastewater treatment facilities. In December 2010, total manageable daily loads (TMDL) were adopted by U.S. EPA in collaboration with the six states and the District of Columbia. These agreements committed to significant reductions of nutrient and sediment loads by 2025, development of locally based watershed implementation plans, and an accountability system including annual milestones and public reporting of progress.
Together, these actions have led to modest progress toward improved water quality and changes in phytoplankton ecology in the bay 16 , although additional nutrient reductions must be reached to decrease APPP and change trophic classification. Our analyses of long-term trends showed flow-adjusted TN and NO 2 + NO 3 loadings doubled from 1945 to 1981, followed by decreases of 19.2% and 5.3% from 1981 to 2012 16 . The slow, upward trajectory of flow-adjusted chl-a for the MH salinity zone is consistent with shallow, downward trends of TN and NO 2 + NO 3 loadings in recent years 16,18 . We point out that simulated 40% reductions of euphotic-layer chl-a or nutrient loadings exceed actual progress since the 1980s, explaining the continuing increases of APPP based on mean-flow model predictions of NPP (Fig. 11a-c).
Decadal contrasts of NPP and APPP in the 1960s and 1970s (Fig. 12a-f ) reflected a combination of water-quality management and climatic effects: (1) lower inputs of bio-optically active constituents in the 1960s accompanied a sequence of low-flow, "dry" years compared to the 1970s, reducing light-limitation for OH, MH, and PH salinity zones and enhancing NPP and APPP; (2) removal of orthophosphate (PO 4 3− ) from detergents enhanced P-limitation in the OH salinity zone, leading to increased N-throughput to MH and PH salinity zones, Figure 12. (a-f) Mean, monthly log 10 NPP (left y-axis) and log 10 chl-a (right y-axis) in the 1960s and 1970s for OH, MH, and PH salinity zones. NPP derived using gam5 ( www.nature.com/scientificreports www.nature.com/scientificreports/ and reductions of NPP and APPP from the 1960s to the 1970s; (3) model predictions of NPP for low-flow, "dry", mean-flow, and high-flow, "wet" conditions were based on predictor variables for flow and salinity that adjusted for climatic effects, with mean-flow predictions of NPP and APPP reflecting these adjustments in the 1960s and 1970s.
Model predictions of NPP from historical reconstructions and recent years led to comparable estimates of APPP for the MH salinity zone. Mean-flow model predictions of NPP produced estimates of APPP >500 g C m −2 y −1 from 2010 to 2015 (Fig. 11b), similar to APPP for the same salinity zone in the 1960s (Fig. 12b). Analogous estimates for OH and PH salinity zones showed a similar pattern, with APPP in the 1960s higher than in recent years. APPP was lower for all three salinity zones in the 1970s than in the 1960s (Fig. 12d-f). These observations and predictions provide historical context for comparison with contemporary conditions, suggesting APPP today is not appreciably different from past rates. We found evidence of lower APPP for MH and PH salinity zones in    Table 4. Statistics for generalized additive models (GAM) of net and gross primary production (NPP, GPP) in Chesapeake Bay using predictor variables (Table 3) as detailed for gam2 (Table 5). GCV = generalized crossvalidation score; AIC = Akaike information criterion; RMSE = root mean square error. www.nature.com/scientificreports www.nature.com/scientificreports/ contemporary estimates, and sensitivity of APPP to reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings. These findings show promise for future reductions of APPP in response to improvements of water quality that would be required to change trophic classification.

Summary
Statistical models of NPP and GPP were developed for Chesapeake Bay with adjustments for climatic effects on hydrology, calibrated with 100 s of shipboard measurements from 1982 to 2004. Model predictions of NPP based on water-quality properties in the 1960s and 1970s and 1985 to 2015 as inputs supported computations of APPP. Simulated reductions of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings led to decreased APPP sufficient to change the trophic classification of the bay. Summarizing: • Statistical models of NPP and GPP calibrated with long-term data included explicit terms to adjust for climatic effects. • These models supported predictions of NPP using historical and monitoring data as predictor variables.
• Model predictions of NPP using historical (1960s, 1970s) and monitoring (1985 to 2015) data as predictor variables supported computations of APPP. • Simulated 40% decreases of euphotic-layer chl-a or TN and NO 2 + NO 3 loadings reduced APPP and changed trophic classification. • Improved water quality is attainable with a reversal of nutrient over-enrichment of the bay sufficient to decrease phytoplankton biomass, but progress to date has been modest compared to goals, exemplified by continuing, high APPP in recent years.  Table 5. GAM models of net primary production (NPP) for the oligohaline (OH) salinity zone. a gam models of NPP for MH and PH salinity zones had the same structures as these models, using input data PP_MH and PP_ PH; models based on the all calibration data (1982 to 2004) used input data PP_ALL and added a categorical term for Salzone (OH, MH, PH); b gam models to predict GPP had the same structure as those for NPP with log_Pbopt_gross substituted for log_Pbopt_net; b gam2 models applied to water-quality data from 1985 to 2015 were used to predict NPP from data files WQ_OH, WQ_MH, and WQ_PH; the explanatory variable Seq_year was omitted from these models to avoid extending trends in calibration data (1982 to 2004) to years outside that time frame; c Substituted log_Euchl with log_Chl in gam models to predict NPP to test models using input data that lacked euphotic-layer chl-a; d Annual TN and NO 2 + NO 3 loadings were replaced with monthly loadings in gam models to predict NPP for historical data, based on the lack of data on Kpar or monthly loadings for the 1960s and 1970s; e Season was omitted as a categorical variable in gam models to predict NPP due to small sample sizes for historical data.   Table 6. Statistical properties and predicted net primary production (NPP) for water-quality data from 1985 to 2015. NPP FIT = model-fitted NPP using observed flow and salinity; NPP MNS = flow-adjusted NPP using long-term means of freshwater flow and salinity.