Microplankton Interactions with Decadal-Scale Nutrient Enrichment in a Deep Estuary, with Implications for Eutrophication—Related Ecosystem Stressors

Long-term (decadal) records of microplankton provide insights into how lower trophic levels of coastal ecosystems respond to nutrient enrichment, over and above shorter-term variability. We used a 15-year seasonal census in the Firth of Thames, a deep, nutrient-enriched estuary in northeastern Aotearoa/New Zealand, to determine microplankton responses to enrichment. Kendall trend analyses showed that dissolved inorganic nitrogen (DIN) and dissolved organic nitrogen were enriched by 99% and 34%, respectively, over 15 years, while phosphorus changed little. Larger phytoplankton (> 2 µm) increased by 46%, including 57% increases by diatoms (mainly large centrics with 93% increase) and nanoflagellates (151% increase). Dinoflagellates decreased by 46%, such that the community shifted from dinoflagellate to diatom/nanoflagellate dominance. Within phytoplankton ≤ 2 µm, picoprokaryotes increased by 369%, while picoeukaryotes changed little. Among microheterotrophs (< 200 µm), bacteria increased by 89%, and small oligotrichs increased by 53%. Trend analyses and multivariate general additive modelling showed that microplankton biomass responded primarily to increased DIN over the 15-year period and secondarily to stratification variation at shorter time scales. The changed biomasses and community composition are explained as responses to increased N:P and food-web interactions. Deleterious changes included increased toxic Pseudo-nitzschia abundance and potentially reduced nutritional quality of the phytoplankton community for grazers. The increased N and larger diatoms indicated potential for increased deposition to sediments, possibly explaining previous observations of lowered denitrification in the Firth during the time series period. The results indicated a continuation of enrichment the Firth has received over decades, with implications for expression of ecosystem stressors of acidification and hypoxia.


Introduction
Coastal ecosystems are threatened by stressors arising from nutrient enrichment and eutrophication (Alexander et al. 2017;Rabalais et al. 2009), including depleted oxygen (O 2 ), lowered pH (Conley et al. 2009;Gobler et al. 2014;Paerl 2006), and changed phytoplankton biodiversity including increases in harmful algae (Capriulo et al. 2002;Sutula 2011). While it is well-known that auto-and heterotrophic microplankton (< 200 µm) are pivotal in mediating water column biogeochemical cycles (Falkowski et al. 1992), less well-understood is how changing microplankton biomass, taxonomic structure, and diversity interacts in expression of these biogeochemical stressors (NRC 2000). When conducted in tandem with biogeochemical studies, long-term (decadal-scale) studies of marine plankton communities can provide useful insights into how lower trophic levels of coastal ecosystems respond to nutrient enrichment and are involved in the expression of eutrophication-related stressors (Eriksen et al. 2019;Kemp et al. 2005;NRC 2000;Philippart et al. 2000).
Communicated by James L. Pinckney The degree to which primary production is exported and mineralised, as well as its cycling within the microbial loop, affects its role in mediating eutrophication-related stressors (Capriulo et al. 2002). Nutrient enrichment can be expected to promote larger-celled phytoplankton, by relieving nutrient limitation in large cells that have small surface-to-volume ratios (Chisholm 1992). Cooper (1995) and Ragueneau et al. (1994) described increases in ratios of large centric diatoms relative to smaller pennate diatoms during early phases of eutrophication of Chesapeake Bay (Kemp et al. 2005). Philippart et al. (2000) showed that enrichmentdriven increases in large cells can lead to accentuated vertical export upon seasonal nutrient limitation and senescence, when cells settle to depth and are mineralised. Enrichment effects on taxonomic structure thus have implications for vertical biomass distribution and export toward mineralisation pathways. Also, as biomass increases due to enrichment, so may the incidence of harmful algal bloom species (Anderson et al. 2008;Cloern 2001) with the potential for associated deleterious ecological effects.
Microheterotrophic (< 200 µm) communities (bacteria, small heterotrophic flagellates and microzooplankton) respond quickly to enrichment of autotrophic populations (Azam et al. 1983;Krstulović et al. 1998). Bacteria drive mineralisation pathways and production of dissolved organic matter originating in large part from phytoplankton, which in turn can support small heterotrophic and mixotrophic flagellates (Antia et al. 1991;Flynn et al. 2019). Microzooplankton grazers adapt quickly to changing bacterial and phytoplankton populations, with the sizes of microzooplankton grazers directly linked to the size of their prey (Hansen et al. 1994). We would therefore expect changes in enrichment and autotrophic population structure to be reflected in microheterotrophs, which will contribute to the respiratory metabolism of the water column, and hence potential for reduced O 2 , increased CO 2 and acidification (Shen et al. 2019;Sunda and Cai 2012;Wallace et al. 2014;Zhou et al. 2020). Concomitant assessment of both autotrophic and heterotrophic components will thus improve understanding of ecosystem responses to enrichment at the base of the food web, with implications for potential ecological stressors.
In this study, we describe decadal-scale changes in the autotrophic and heterotrophic plankton in the Firth of Thames, in northern Aotearoa/New Zealand (NZ), over 15 years (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). The Firth is a large (1100 km 2 , 15 m average depth), drowned-valley estuary (Hume et al. 2016) that opens to the Hauraki Gulf (Fig. 1). It supports some of NZ's largest marine farms (for mussels) and is an important region for recruitment to its largest inshore fin-fishery (for snapper) (Kelly et al. 2020;Peart 2016;Zeldis and Francis 1998). Total nitrogen (TN) loading to the Firth, primarily from the Waihou and Piako rivers at its head (Fig. 1B), has increased substantially above natural levels (by ~ 82%: Snelder et al. 2018). The agricultural intensification driving this loading increase is historical, having occurred over many decades, commencing in the early 1900s with drainage of wetlands and peat swamps, but undergoing rapid expansion in the 1980s and 1990s primarily for dairy farming (Judd 2015;Kelly et al. 2020;Peart 2016). Currently, about 78% of TN runoff to the Firth is anthropogenic (Vant 2013), and catchment runoff now dominates its nutrient loading, with ~ 85% of its dissolved inorganic nitrogen (DIN) supply arising from its catchment (Zeldis and Swaney 2018). As a consequence, the Firth now sustains volumetric primary production roughly double that of the seaward Hauraki Gulf (Gall and Zeldis 2011;Zeldis et al. 2022;Zeldis and Willis 2015). It has a pronounced seasonal cycle of net ecosystem metabolism that evolves from net autotrophy in spring and early summer to net heterotrophy from late summer through autumn, when it often expresses significantly depressed pH and hypoxia (Law et al. 2019;Zeldis et al. 2022;Zeldis and Swaney 2018).
Nutrient enrichment of the Firth has continued during the more recent period of the present study (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013), with increases in DIN and dissolved organic nitrogen (DON) concentrations (Zeldis and Swaney 2018). To investigate how microplankton at the base of the Firth food web have responded to nutrient enrichment, here we describe decadalscale changes in its phytoplankton and microheterotrophic biomasses and diversity. We place these changes into ecological context by assessing the inter-annual and seasonal dynamics of the microplankton taxonomic groups, their responses to the changed nutrient environment and trophic interactions. We assess potential causes of the N increases the Firth has received over the recent period and the role of microplankton in that enrichment. We conclude by describing the roles of microplankton in the expression of pH and O 2 -related stressors currently seen in the Firth.

Sampling of Nutrients, Physical Properties, Phytoplankton, Bacteria and Microzooplankton
Surveys were conducted at the NIWA (National Institute of Water and Atmospheric Research, New Zealand) Firth of Thames monitoring site (Fig. 1B) during 57 voyages made at seasonal (quarterly) intervals from September 1998 to December 2013. Sampling was typically in austral spring (September/October), summer (December), autumn (March) and winter (July), except for July 2001 to December 2002, when sampling was not conducted. Physical samples were acquired from September 1998 to December 2013 using a conductivity-temperature-depth (CTD: Seabird 911plus) instrument, fitted with O 2 and water transmissivity sensors and a water sampling rosette fitted with 10-L Niskin bottles, deployed from the RVs Kaharoa (28 m), Tangaroa (78 m) or Star Keys (19 m). Physical samples were resolved at 1m depth intervals through the 40 m water column. CTD sensors were factory calibrated within the interval recommended by the manufacturer (2 years) and had accuracies of 0.002 °C, 0.003 mS/cm and 2% of surface O 2 saturation, respectively. Transmissivity was sampled using a Seabird C-Star transmissometer (25 cm path length).
Additional seasonal surveys were conducted on a station transect through the Firth and into the Hauraki Gulf ( Fig. 1B) in spring, summer, autumn and winter 2012-2013 and autumn 2010, as part of research in the wider Hauraki Gulf/Firth region that examined its carbonate and O 2 environments (Zeldis et al. 2022). These surveys used CTDrosette sampling and analysis protocols identical to those at the Firth monitoring site.
Samples (250 mL) for phytoplankton (> 2 µm in size) and microzooplankton (< 200 µm) were acquired from September 1998 to July 2013. Samples were preserved in Lugol's iodine (1% final concentration) on board the vessel. In the lab, 50-250 mL sub-samples were settled in measuring cylinders with settling times of 0.5 h per millimetre of settling column depth. The supernatant was then siphoned off, and the remaining 10 mL was transferred to Utermohl chambers and resettled. Phytoplankton > 2 µm were then enumerated using a Leica DMI3000B inverted microscope at × 100 to × 600 magnification (Edler and Elbrächter 2010). The cell-size lower detection limit for these microscope counts was ~ 2 µm. Samples were counted as soon as practicable after collection. Although samples were occasionally held for longer durations (up to months), there was no trend in how long they were held before analysis, within the ~ 15-year time series. All abundant phytoplankton > 2 µm at concentrations > 100 cells per L were identified to genus or species level where possible, counted and assigned to algal groups including diatoms (centric and pennate) and dinoflagellates (naked and armoured). Known toxin-producing species of dinoflagellates and diatoms were itemised into additional subgroups. Small flagellated cells (excluding dinoflagellates) were enumerated as nanoflagellates (most < 20 µm), including unspeciated, Prymnesiophytes, Cryptophytes, Chrysophytes, unidentified monads and rarer populations of identified silicoflagellates, Raphidophytes and euglenoids. Heterotrophic and mixotrophic nanoflagellates were not separated from autotrophic nanoflagellates, and their populations were assessed as an undifferentiated group. All phytoplankton > 2 µm sample analyses were conducted by the first author. Any changes in nomenclature were retrospectively applied to older data to keep systematics consistent. Acidified Lugol's iodine used for preservation of phytoplankton > 2 µm was likely to have eliminated any calcium-based algae (coccoliths). However, subsequent counts with alkaline Lugol's iodine showed this group was not an important component of biomass in these waters.
Microzooplankton, including ciliates and larval copepods (nauplii), were determined from the same bulk Lugol's iodine samples as phytoplankton > 2 µm. Samples were analysed following similar protocols to phytoplankton, but settled volumes were 250 mL with supernatant siphoned and the remaining ~ 20 mL then transferred to Utermöhl chambers and resettled (Edler and Elbrächter 2010). A wild inverted microscope was used for enumeration during 1998-2001 and then subsequently a Leitz Fluovert FS inverted microscope.
Microzooplankton were identified to genus where possible but with no differentiation of plastidic ciliates. All microzooplankton analyses were conducted by the fourth author.
Picophytoplankton (phytoplankton ≤ 2 µm) and bacteria were sampled between December 2002 and December 2013 and frozen immediately on board in liquid N 2 in 1.8-mL cryovials (in triplicate). In the lab, samples were stored at −80 °C and were thawed immediately before counting on a FACSCalibur 15 mW air-cooled argon-ion laser flow cytometer (fixed at 488 nm: Becton Dickinson, Mountain View, CA, USA) with CellQuest Vers. 3.3 software. Counts for picoprokaryotes (Synechococcus), picoeukaryotes and bacteria (stained) followed the methods of Hall et al. (2006). All flow cytometric samples were analysed by the sixth author within 2-3 months of collection.
Phytoplankton > 2 µm biovolume estimates were calculated from the dimensions of each taxon using approximated geometric shapes (spheres, cones, ellipsoids), initially following Hillebrand et al. (1999) with some modifications following Sun and Liu (2003) and Olenina et al. (2006). Any changes in shape estimates for taxa were applied retrospectively, to avoid bias in estimates with time. Size categories within taxa were also used to reduce bias (e.g. 'small Gyrodinium' and 'large Gyrodinium'). The biovolumes were used to calculate cell carbon (mg C m −3 ) using equations from Montagnes and Franklin (2001) for diatoms and Strathman (1967), Menden-Deuer and Lessard (2000) and Putt and Stoecker (1989) for dinoflagellates and nanoflagellates.
Ciliate cell C was estimated from dimensions of 10-20 randomly chosen individuals of each taxon using approximate geometric shapes and then converted to C biomass using a factor of 0.19 pg C µm −3 (Putt and Stoecker 1989). Copepod naupliar cell C was estimated using a conversion factor of 0.032 µg C individual −1 (Uye and Sano 1998). Cell C biomasses for picoeukaryotes (1319 fg C cell −1 ) and picoprokaryotes (154 fg C cell −1 ) were estimated following Buitenhuis et al. (2012). Bacteria C biomass was estimated following Fukuda et al. (1998) using 30.2 fg C cell −1 .

Data Analysis
Data for physical parameters, nutrients and biological groups were depth-integrated using: where I denotes a sample level number (incrementing from 1 at the sea surface), B i denotes physical values or concentrations (depending on variate type) at sample level I and ΔZ i,i+1 is the thickness (m) of the layer between sequential sampling depths. In analyses where data were grouped by seasons, seasons were defined as September to November, December to February, March to May and June to August, for austral spring, summer, autumn and winter, respectively. Interpolated displays of data were visualised using Ocean Data View software (Schlitzer 2013).
Time-series temporal trends were assessed using nonparametric seasonal Kendall trend tests using 'Time Trends' software (Jowett 2019) with trend slope directions tested at p = 0.05. Interpretation of the trends at that level considers the following precepts from McBride (2019) and Jowett (2019): (1) if the sign of both the upper and lower 90% (onesided) confidence limits for the time trend (Sen) slope are the same (or one limit is zero), then it can be inferred that there is a trend at Kendall p = 0.05, and (2) while the Kendall p value for significance of slope indicates the strength of trends, p values ≥ 0.05 are not sufficient evidence to suggest that no trend exists, but only that the null hypothesis (i.e. there is no trend) is 'not rejected' (McBride 2019). To accommodate this, and following McBride (2019) and Jowett (2019), the probability that the slope is truly below or above zero is calculated and further assessed for likelihood using the terms in Table 1 (Jowett 2019;Mastrandrea et al. 2010). Full trend statistics are given in the Supplementary Material Tables S1-S4 of this paper.
To assess changes of nutrient and microplankton values from the start to end of their respective time series, their median values (V) at the start and end were evaluated using: where V M is the median value at the time series mid-point, S is the Sen slope of the Kendall seasonal trend standardised to a year and expressed in the units of the variable (Jowett 2019) and L is the time series length (y). Percentage changes from start to end of the respective time series were expressed as (V end -V start )/V start *100 and presented for values that had 'increasing/decreasing trend possible' or stronger. Time series lengths were taken as 15.5 years for nutrients and pigments, 15.0 years for phytoplankton > 2 µm and microzooplankton and 11.25 years for picophytoplankton and bacteria, based on their respective start and end dates ('Sampling of Nutrients, Physical Properties, Phytoplankton, Bacteria and Microzooplankton' section). Associations between depth-integrated physical and nutrient predictors and phytoplankton and bacteria biomasses (i.e. the broad taxonomic groups: diatoms, dinoflagellates, nanoflagellates, picoprokaryotes, picoeukaryotes, oligotrichs, tintinnids, nauplii and bacteria) were assessed with multivariate General Additive Models (GAMs), using the package 'mgcv' (Wood 2006(Wood , 2011. GAMs were used to assess the potential for linear and non-linear covariance between continuous physical (temperature, salinity, stratification) and nutrient (DIN, DON, DIP) predictors and the categorical predictors 'group' (diatoms, dinoflagellates, nanoflagellates, picoprokaryotes, picoeukaryotes, oligotrichs, tintinnids, nauplii and bacteria) and 'season' (spring, summer, autumn, winter). Assumptions of normality (Q-Q plot) and homogeneity of variance (Levene's test), as well as 'concurvity' for general additive model analysis (an estimate of redundancy among explanatory variables), were checked for the multivariate model. GAM models were fitted with 'cr' (cubic regression) splines, using a k-value of 10 (i.e. the number of 'knots' denoting the complexity of the non-linear fit), and the distribution family 'tweedie'. Selection procedures were implemented to penalise and remove factors with poor explanatory power.

Dissolved Nutrient Trends
Concentrations of DIN plotted by depth and time over the 15-year time series at the Firth monitoring site showed that the upper water column (< 20 m) sustained persistently depleted DIN median concentrations (~ 0.5 mmol m −3 ), especially from spring through autumn, with higher concentrations below 20 m (Fig. 2). These DIN differences with respect to depth corresponded with observations by Zeldis and Swaney (2018) who showed (their Fig. 9) that the Firth pycnocline typically formed around 20 m depth, about mid-depth in the 40 m water column, separating the water column into a deeper layer (lower, ≥ 20 m) isolated from the upper layer (upper, < 20 m) by density discontinuity. DIP also was often depleted in the upper water column ( Fig. 2) but less so than DIN. Overall, DIN:DIP was low, with median values < 3.5, well below the typical 16:1 Redfield ratio for N and P balanced growth (Redfield et al. 1963), indicating strong N-limitation of phytoplankton growth.
Seasonal Kendall trend tests (Fig. 3) showed that water column-integrated DIN had a 'very likely' increase of 4.3% year −1 (Table S1A), with increases of 4.1% year −1 and 4.4% year −1 in the lower and upper water columns, respectively ('likely' and 'very likely' trends) (Table S1B). Analysis by season showed that the season with the highest likelihood of an increasing DIN trend was summer (6.0% year −1 : Fig. 3; Table S1C). The NH 4 + -N component of DIN had 'likely increasing' trend in the upper water column. DON was abundant relative to DIN (Table S1A) and had 'very likely' increase at 1.9% year −1 , with the highest likelihood of increasing trend in autumn (2.3% year −1 ). In contrast, whole column-integrated DIP had no important change and a 'likely' decrease (− 2.1% year −1 ) in the upper column. DOP was less abundant than DIP (Table S1A) and had a Overall, during their 15.5-year time series, whole water column-integrated DIN and DON had increases of 99% and 34%, respectively. DRP and DRSi changed little and DOP decreased by 54%.

Autotrophic Trends
The biomasses of phytoplankton > 2 µm and picophytoplankton (< 2 µm) groups over their respective time series at the monitoring site were generally the highest in the upper (< 20 m) water column (Fig. 2) that corresponded with persistently depleted DIN concentrations at those depths. On occasion, biomass was distributed more evenly throughout the water column or concentrated in the lower column (≥ 20 m), particularly for diatoms. There was a pattern for dinoflagellates to be more abundant than diatoms early in the time series period, but this pattern reversed later in the period (Fig. 2).
Separate trend analyses of biomasses integrated over upper (< 20 m) and lower (≥ 20 m) water column strata ( Fig. 4A; Table S2B) showed that total phytoplankton > 2 µm biomass increases were stronger in the lower column than the upper column ('very likely' vs 'possible' increasing trends, respectively). Total diatom biomass had a weak 'possible' increasing trend in the lower water column (3.8% year −1 ) and no likely trend in the upper water column. Large centric diatoms, however, had a stronger 'increasing trend likely' (5.1% year −1 ) in the lower column. Pennate diatoms had an 'increasing trend possible' in the lower column and upper column, increasing 2.5% year −1 and 2.9% year −1 , respectively. Toxic diatoms in the lower column had a 'virtually certain' increasing trend (14.5% year −1 ) and a 'very likely' increasing trend (11.2% year −1 ) in the upper column. Total dinoflagellates had no likely trend in the lower column and a 'possible' decreasing trend in the upper column,

Fig. 3
Likelihoods and strengths of temporal trends for nutrient concentrations (mmol m −2 ) determined using seasonal Kendall trend tests. The first arrow in each row shows trend direction and likelihood for concentrations integrated over the 40 m water column using the arrow symbols described in Table 1. Subsequent arrows in each row show trend results for concentrations integrated over upper (< 20 m) and lower (≥ 20 m) water column depth strata and seasons, for which trends were 'possible' or stronger in likelihood for 40 m integrated concentrations. Annotated percentages are percent annual changes of concentrations, and seasons are annotated as spring (Sp), summer (Su), autumn (Au) and winter (Wi) while armoured dinoflagellates had no likely trends. Naked and toxic dinoflagellates had no trends in the lower column but 'virtually certain' and 'very likely' decreases (−6.1% year −1 , −9.5% year −1 , respectively) in the upper column. Nanoflagellates had a 'virtually certain' increasing trend in the lower column (6.7% year −1 ) and 'very likely' increase in the upper column (6.5% year −1 ). Picoprokaryotes had 'virtually certain' increasing trends in both lower (9.8% year −1 ) and upper water columns (11.8% year −1 ), while picoeukaryotes had a 'possible' increase (3.0%) in the lower column.
Median values of the taxonomic group biomasses at their time series midpoints (Fig. 5) showed that diatoms, Fig. 4 Likelihoods and strengths of temporal trends for concentrations of A phytoplankton (mg C m −2 ) and B chlorophyll a and phaeopigment-a (mg pigment m −2 ), determined using seasonal Kendall trend tests. The first arrow in each row shows direction and likelihood of change in total biomass of the groups integrated over the 40 m water column, using the arrow symbols described in Table 1. Subse-quent entries in each row show trend results for biomass integrated over lower (≥ 20 m) and upper (< 20 m) water column depth strata and seasons for which trends were 'possible' or stronger in likelihood, for 40 m integrated biomass. Annotated percentages are percent annual changes of biomass, and seasons are annotated as spring (Sp), summer (Su), autumn (Au) and winter (Wi) dinoflagellates and nanoflagellates contributed roughly equal proportions to total phytoplankton > 2 µm. The diatoms were composed mostly of centric types, dominated by several large genera, including Coscinodiscus, Trieres (formerly Odontella), Rhizosolenia, Proboscia, Thalassiosira and Chaetoceros. Pennate diatoms were less abundant than the centrics and were mainly Pseudonitzschia, Nitzschia and Navicula. The toxic diatoms were Pseudo-nitzschia and were rare overall but occasionally abundant in spring. Dinoflagellates were mostly armoured, followed by naked types, the latter contributing most to the overall declining dinoflagellate biomass trends (Fig. 4A). The declining toxic dinoflagellate group included the genera Karenia, Alexandrium, Dinophysis, Gonyaulax, Lingulodinium and the less common Ostreopsis and cf. Gambierdiscus. Nanoflagellates contributed strongly to the increasing phytoplankton > 2 µm biomass (Fig. 4A) and were dominated by monads, Cryptophyceae, Chrysophyceae and Prymnesiophyceae with occasional euglenoids, silicoflagellates, Octactis/Dictyocha and toxic raphidophytes Heterosigma, Fibrocapsa, Chattonella and Prymnesium. There were only three instances when any of the toxic components of the nanoflagellates contributed more than 0.1% of biomass (not shown in Fig. 5). The picoprokaryotes, while increasing strongly over the time series, contributed relatively little to overall phytoplankton biomass (~ 6% on average) (Fig. 5). In contrast, the picoeukaryotes, while having few important biomass trends (Fig. 4A), contributed substantial biomass, comparable to the summed biomasses of diatoms, dinoflagellates and nanoflagellates.
Chl-a pigment (Fig. 4B), integrated over the whole 40 m water column, had a 'possible' increase (1.7% year −1 ), while its degradation product, phaeo-a, had no important trend ( Fig. 4B; Table S3A). However, chl-a and phaeo-a had 'virtually certain' and 'very likely' increases in the lower water column, (4.9 and 3.1% year −1 , respectively). In the upper column, chl-a had no trend, and phaeo-a had a 'likely' decrease. Seasonal patterns for pigments integrated over the whole water column were weak with only 'possible' increases in winter for chl-a and in autumn for phaeo-a.
The increases of pigments in the lower water column reflected the increases seen in total phytoplankton > 2 µm, particularly centric diatoms and nanoflagellates (Fig. 4A). While chl-a integrated over the whole water column increased by only 31% and phaeo-a had little change, chla and phaeo-a integrated over the lower water column (≥ 20 m) increased by 123 and 63%, respectively.

Drivers of Biomass
Multivariate GAM analysis identified DIN concentration and stratification as having the strongest relationships with phytoplankton and bacterial biomasses, with DIN identified as the dominant driver and stratification secondary as shown by the higher F ratios and p values for DIN (Table 2). Concurvity or collinearity of variables showed that model parameters were between 0.72 and 0.76 (values close to one indicate redundancy), with DIN showing the least covariance with other parameters (Table 2). Bacteria, diatoms, picoeukaryotes and nanoflagellates all had increasing linear or saturating relationships with increasing DIN (Fig. 7), while dinoflagellates and picoprokaryotes had decreasing relationships. Factors with poor explanatory power that were dropped from the analysis were DON, small and large oligotrichs and nauplii. Lagging the GAM model by one season (i.e., ca 3 months) with respect to biomass (not shown) yielded poorer fits than the non-lagged model, with adjusted r 2 dropping from 0.63 to < 0.10. Given the likely quick responses of biomass to changes in nutrient and physical conditions (days rather than months), this is not unexpected.

Fig. 6
Likelihoods and strengths of temporal trends for heterotrophs (mg C m −2 ), determined using seasonal Kendall trend tests. The first arrow in each row shows direction and likelihood for biomass of the groups integrated over the 40 m water column, using the arrow symbols described in Table 1. Subsequent arrows in each row show trend results integrated over lower (≥ 20 m) and upper (< 20 m) water column depth strata and seasons for which trends were 'possible' or stronger in likelihood, for 40 m integrated biomasses. Annotated percentages are percent annual changes of concentrations, and seasons are annotated as spring (Sp), summer (Su), autumn (Au) and winter (Wi) Table 2 General Additive Model (GAM) output for assessing the influence of water column-mean values of physical properties, temperature (°C), salinity (psu), density stratification (Brunt-Väisälä frequency: cycles h −1 ) and water column-integrated nutrient concentrations (mmol m −2 ) on microbial plankton biomasses (mg C m −2 ) over seasonal surveys 1998-2013. Parameters significantly related to biomass (p < 0.05) are highlighted in bold. The table includes diagnostics of the relationship between each parameter and biomass (Edf estimated degrees of freedom, Ref.df reference degrees of freedom, F and p value) and the concurvity estimate for each parameter (Wood 2006). The overall model fit had an adjusted r 2 of 0.63 and 62.4% deviance explained. Categorical comparisons between categorical factors 'microbial plankton group' and 'season' are given in The interactive effects of DIN and stratification on biomass ( Fig. 8) showed that under strong (high values) water column stratification and low DIN concentrations, biomass was highly suppressed. However, increases in DIN were associated with strong increases in biomass across the full range of stratification values. Overall, the relationship predicted a saturation of biomass increase at DIN concentrations between 100 and 150 mmol m −2 , although inferred declines in biomass beyond these concentrations were described by increasingly sparse data points (Fig. 8). To further isolate the relative effects of DIN and stratification on biomass, inter-annual and intra-annual (seasonal) variations in these parameters were considered. While DIN was shown to have a 'very likely' increasing inter-annual trend (Fig. 3, Table S1), water column temperature, salinity and stratification all had 'unlikely' trends (Table 3). However, stratification varied significantly intra-annually (between seasons), while DIN had no significant seasonal variation (Table 4). These results indicated that the DIN effect on biomass detected by the GAM arose primarily via its strong increasing trend over the 15-year time series, whereas the relatively weak effect of stratification on biomass arose at shorter time scales, potentially via stochastic, year-to-year variability and intra-annual seasonal variation.
Comparison of categorical factors 'microplankton group' and 'season' used in the GAM to predict effects of nutrient and physical variables on biomass showed that effects due to factor 'microplankton group' were stronger (p < 0.001) than effects due to factor 'season', which were largely nonsignificant (Table S5). The relatively strong effects of 'microplankton group' likely arose from the contrasting responses of diatom, nanoflagellate, picoeukaryote and bacteria groups (increasing biomass) and the dinoflagellate group (decreasing biomass) to increasing DIN (Fig. 7). The relatively weak response to factor 'season' was likely due to lack of biomass variations across seasons which were largely non-significant, apart from that of picoeukaryotes (Table 4).
The seasonal surveys made in 2012-2013 and in 2010 which traversed the Firth (Fig. 9) showed that chl-a and phaeo-a decreased in the upper column (< 20 m) and Fig. 7 Relationship between DIN concentration (mmol m −2 ) and water column-integrated biomass (mg C m −2 ) of phytoplankton taxonomic groups and bacteria. Regressions are fitted using General Additive Models (± 0.95 CI). x-and y-axis are square root-transformed increased in the lower column between spring and summer. These patterns were reflected by increased transmittance in the upper water column and decreased transmittance in the lower column. Nitrate was depleted to low levels (≤ 0.5 mmol m −3 ) in the upper water column from spring through autumn, accompanied by strengthening stratification spring to summer. The deeply distributed biomass in summer was accompanied by increased NH 4 + -N and decreased O 2 , indicating that it was senescing and mineralising. Interannual variation in the seasonal occurrence of these processes was seen by comparing the surveys in summer 2013 and autumn 2010 (Fig. 9) that showed similar chl-a, phaeoa, transmittance, stratification and O 2 patterns but with a seasonal shift.

Inter-annual and Seasonal Dynamics of Microplankton Biomass
The 15-year trend analyses showed increasing DIN coincided with increases in several auto-and heterotrophic group biomasses, while the GAM analysis identified DIN increase as the dominant correlate of biomass changes, with stratification secondarily important. While having a strong increasing trend over the time series, DIN trends did not vary significantly between seasons (Table 4), indicating that it affected biomass primarily via its longterm trend over the time series, across all degrees of Fig. 8 Interactive response of phytoplankton/bacteria biomass (as predicted by General Additive Models: green = low, orange = high biomass) to DIN concentration (mmol m −2 ) and water column stratification (Brunt-Väisälä frequency: (cycles h −1 : higher values mean more strongly stratified)). Fitted mesh is constrained to interpolate 25% beyond actual data Table 3 Seasonal Kendall time trend tests for water column median temperature (°C), salinity (psu) and density stratification (Brunt-Väisälä frequency: cycles h −1 ). Shown are numbers of seasons sampled, the median value over 1998-2013 time series, the significance level of trend slope (Kendall P), the % annual change, the Sen slope, 90% confidence limits of the slope, the probability that the slope was truly above or below zero and confidence of trend direction  (Fig. 8). This association with DIN was consistent with findings in coastal ecosystems showing that biomass is demonstrably related to nutrient loading rates and concentrations (Boynton and Kemp 2008;Boynton et al. 1982;Monbet 1992), albeit with variability in the relationships across different types of coastal ecosystems (Cloern and Jassby 2008). The long residence times and relatively clear waters of large coastal embayments such as the Firth mean such systems are sensitive to nutrient enrichment (Zeldis et al. 2022). They are susceptible because, in addition to having potentially high nutrient loading from large catchments, they can support phytoplankton blooms which undergo complete cycles of growth, consumption and senescence, prior to seaward export (Ferreira et al. 2005). Stratification had a significant but weaker effect than DIN, with increased stratification depressing microplankton biomass. This could be caused by decreased vertical flux of nutrients from the pycnocline into the mixed layer with increased stratification, an effect considered an important controlling factor in primary production in seasonally stratified coastal waters (Eppley et al. 1979;Mann and Lazier 1991). In the Firth, this may have operated at both inter-annual and intra-annual (seasonal) time scales, affecting biomass at both frequencies. Although stratification strength had no consistent trend over the 15-year time series (Table 3), it did vary from year to year at the Firth monitoring site, as shown by Zeldis and Swaney (2018) (their Fig. 9). This could have accounted for its effects on biomass  in the GAM, albeit without a consistent trend over the whole time series. Stratification also varied strongly between seasons at the monitoring site (Table 4), as was also shown by seasonal surveys made along the length of the Firth (Fig. 9).

Linking Changes in the Microplankton Community to Nutrient Ratios
Changes in autotrophic and heterotrophic groups over their time series were accompanied by increasing trends in concentrations of DIN and DON at the Firth monitoring site. In contrast, DIP, DOP and DRSi had no trends or decreasing trends, resulting in changed ratios among the nutrient species. Nutrient ratios can affect phytoplankton composition because different taxonomic groups display differential preferences for specific nutrient ratios or nutrient forms, as proposed in the 'nutrient ratio' hypothesis (Smayda 1990(Smayda , 1997Tilman 1977). This suggested that the increase in the ratios of N to P and DRSi may have underpinned changes in biomass ratios among taxonomic groups. Increased diatom biomass and decreased dinoflagellate biomass during 1998-2013 (Figs. 2 and 4A) altered the ratio between these groups across the time series. The dominance of dinoflagellates at the start of the time series in 1998 resembled that sampled 2 years prior, in 1996-1997 (Chang et al. 2003), who reported a dominance of dinoflagellates over diatoms in the inner Hauraki Gulf. Dinoflagellates have a greater P requirement to support their maximum photosynthetic and respiration rates than other phytoplankton groups (Reynolds 2006). Increased N and stable or declining P resulted in an increasing N:P ratio that could have favoured other, less P-dependant groups, including diatoms. DOP concentrations declined, which also may have been a limiting factor because dinoflagellates can rely on DOP for nutrition (Lin et al. 2015). Silica was also sufficient for diatom growth, with Si:N ratio ~ 5:1, being far from limiting for diatoms (~ 1:1: Chang et al. (2003); (Sarthou et al. 2005). The increasing biomass trends for larger diatoms were consistent with observations that N-stress relief during eutrophication can favour the growth of larger cells, based on cellular surfaceto-volume considerations (Chisholm 1992;Philippart et al. 2000;Ragueneau et al. 1994;Stolte and Riegman 1995).
Nanoflagellates and bacteria also increased substantially over the time series, likely due to a combination of factors. Along with increased DIN supply, there was sufficient Si, which would benefit silico-flagellates and small Chrysophytes, which require adequate Si levels. Increasing DON would have also supported groups which can consume dissolved organic matter, like bacteria and some mixotrophic and heterotrophic flagellates, that can take up DON directly or via grazing (Antia et al. 1991;Berman and Bronk 2003;Flynn et al. 2019;Granéli et al. 1999;Sanderson et al. 2008;Seitzinger et al. 2002). The substantial increases of bacteria and nanoflagellates over the time series were strong in summer and coincided with maximum increases in DIN, suggesting a strong trophic response was occurring through these groups.

Trophic Linkages Within the Microplankton Community
As well as microplankton community changes directly related to dissolved nutrient supply, other changes appeared related to trophic food-web linkages. Small heterotrophic and mixotrophic flagellates grazers are known to respond quickly to increases in bacteria and picoplankton (Capriulo et al. 2002;Duarte et al. 2000;Edwards 2019;Stoecker et al. 2017) and to act as primary consumers of those prey (Duarte et al. 2000;Flynn et al. 2019;Hall et al. 1993;Hansen et al. 1994;Safi and Hall 1999). Although we did not separately enumerate autotrophs, heterotrophs and mixotrophs in our nanoflagellates group, we assume that their strongly increasing biomass indicated that heterotrophy and mixotrophy were important in that group, as found in other coastal systems in New Zealand and elsewhere (Flynn et al. 2019;Riemann et al. 1995;Safi and Hall 1999;Unrein et al. 2007). Over the time series, DON increased most strongly in autumn, potentially from mineralisation of organic matter in summer and autumn ( Fig. 3; Zeldis et al. (2022)). This mineralisation would be driven by bacteria, which increased strongly in summer, over the time series.
Naked dinoflagellates can also be important mixotrophic grazers (Lin et al. 2015) but their numbers declined. This may have been linked to their generally larger size. Mixotrophic dinoflagellates often have a 1:1 ratio in size, with respect to their prey (Hansen et al. 1994), indicating they would predominately consume mid-size prey (~ 20-40 µm). Prey of that size (i.e., Phytoplankton, ~ 20-40 µm) did not increase over the time series, unlike smaller nanoflagellates (dominated by 4-10 µm cells), bacteria, prokaryotic picoplankton (< 2 µm) and oligotrichs < 20 µm in size. Increased biomass of small (< 20 µm) oligotrichs, in contrast, was potentially supported by the increased bacteria and picoprokaryotes (Christaki et al. 1999). Small oligotrichs graze primarily upon these very small prey, unlike larger (> 20 µm) oligotrichs and tintinnids, which had little change or declined over the study period. Ciliates have an optimal prey-to-predator size ratio of about 18:1 (Hansen et al. 1994) indicating that while bacteria and picoplankton may have occupied a suitable prey size range for the smaller oligotrichs, they would not have been optimal for larger ciliates.
Interestingly, picoeukaryotes had a positive response to higher DIN concentration in the GAM analysis ( Fig. 7) but little inter-annual trend over the time-series (Fig. 4A). Their response to DIN within the GAM was therefore likely due to their relatively strong intra-annual (seasonal) response (Table 4), with wide variation across the seasons. The lack of inter-annual response by picoeukaryotes may be the result of tighter coupling between this population and nanoflagellate and oligotrich grazers, which often dominate grazing upon picoeukaryotes (Rassoulzadegan et al. 1988;Samuelsson and Andersson 2003). Picoeukaryotes are likely preferred by nanoflagellates and oligotrichs over picoprokaryotes or bacteria, because of their greater nutritional value (Monger and Landry 1991;Worden et al. 2004).
The weak increase of small ciliates and the lack of increase of larger ciliates may also be a function of topdown control by copepods. Strong predation pressure by copepods on ciliates is reported in subtropical waters east of New Zealand Décima 2020), andDolan (1991) reported strong predation pressure on ciliates by copepod nauplii in coastal waters (Chesapeake Bay, USA). These studies indicate that grazing on ciliates may have supported the increasing trend for copepod nauplii found in the present study.

Deleterious Algal Taxa and Phytoplankton Community Nutritional Quality
Nitrogen enrichment and changed nutrient ratios have been linked with increased occurrence of harmful algal blooms (Davidson et al. 2012;Glibert and Burford 2017). In this study, the toxic pennate diatom genus Pseudo-nitzschia spp. increased strongly (~ 4800%), mostly in spring, consistent with Rhodes et al. (2013) who found Pseudo-nitzschia spp. bloomed most commonly in spring and summer in NZ North Island coastal waters. In contrast, toxic dinoflagellate biomass declined, mainly in summer. Consideration of dissolved nutrient ratios may again help to explain these changes. As discussed, dinoflagellate growth is optimised at low N:P ratios (Hodgkiss and Ho 1997), but trends in the Firth indicated increasing N:P ratios. In the Gulf of Mexico, Pseudo-nitzschia spp. cell counts increased in parallel with increased NO 3 − -N loading over several decades (Parsons et al. 2002;Turner and Rabalais 1991). Our survey has shown that Pseudo-nitzschia spp. biomass increased in the Firth, and, if this trend continues, it may become more problematic given the most toxic strain (P. australis) is reported in this region (Rhodes et al. 2013).
The trends of declining dinoflagellates and increasing diatoms which we observed over this time series could indicate a shift in nutritional quality of the phytoplankton for grazers at higher trophic levels. Dinoflagellates and nanoflagellates have been shown to be considerably more nutritious than diatoms (in terms of C and N constituents) on a chl-a or volume-specific basis (Chan 1980;Hitchcock 1982;Rey et al. 2001), suggesting that the higher densities of dinoflagellates at the start of the time series would have been more favourable for secondary production than later. It was previously shown that the planktonic system of the Hauraki Gulf had greater production of mesozooplankton and larval fish in years with higher dinoflagellate than diatom composition (Zeldis et al. 2005). The present results may indicate that N enrichment has set conditions for a less nutritious environment for mesozooplankton in the Firth. This may have been offset, however, by the increases of nanoflagellates and small ciliates, which can be passed quickly through nauplii and larger microzooplankton (Dolan 1991;Samuelsson and Andersson 2003).

Role of Microplankton in Nutrient Enrichment and Eutrophication
The results showed that increased microplankton biomasses and changed taxonomic group composition accompanied the changed nutrient environment (strong increases in N, no change or decreases in P) over the 1998-2013 period. In the case of P, concentrations declined slightly in the Firth river inputs during the period (by ~ 0.6-1.5% year −1 Vant (2013)), likely driven by point source improvements in the catchment, consistent with the flat or decreasing trends in DIP and DOP in the Firth. For N, Zeldis and Swaney (2018) described the balance of offshore-and catchment-derived N supply to the system. They showed that neither offshore physical oceanographic conditions (upwelling) nor enhanced estuarine exchange changed within the Hauraki Gulf/Firth system, indicating these factors were unlikely to be drivers of changed N. Changes in water quality in the Firth river inputs also were not strong drivers of changed N loading, with only small increases in river TN concentrations (maximum 1.7% year −1 : Vant (2013); Zeldis and Swaney (2018)), much less than the strongly increasing trends of DIN and DON over the period. Thus, there was little evidence of either offshore-or catchment-side change during the monitoring period that would explain the N and microplankton biomass increases in the Firth.
Another potential explanation for the observed N enrichment is changed internal nutrient cycling in the Firth sediments, related to denitrification. The nutrient mass-balance analysis of Zeldis and Swaney (2018) found that the denitrification rate in the Firth decreased by ~ 58%, between surveys made in 2000-2001 and 2012-2013. Studies made elsewhere have shown that the denitrification rate can be inhibited by organic enrichment in sediments (Cook et al. 2004;Eyre and Ferguson 2009;Hale et al. 2016;Harris et al. 1996;Kemp et al. 1990;Tait et al. 2020), raising the possibility that organic enrichment of sediments could have decreased denitrification and increased availability of DIN to primary producers, thus driving a positive feedback on eutrophication. This was supported by sediment biogeochemistry studies made at the Firth monitoring site in 2003 and 2012 (Zeldis and Swaney 2018), showing that sediment O 2 demand increased over that period by ~ 1.5 to 1.7-fold, indicating organic enrichment.
Water column N enrichment can enhance downward transport of phytoplankton, leading to enrichment of sediments (Kemp et al. 2005;NRC 2000). This is due not only to higher algal production generally, but also shifts toward larger algal species with higher sedimentation characteristics (Philippart et al. 2000). In the Firth, larger taxa (phytoplankton > 2 µm) were important contributors to total phytoplankton biomass (Fig. 5) and had increasing trends (Fig. 4A). This was particularly the case for large centric diatoms, which had their largest increases in autumn, toward the end of the production season. Nanoflagellate biomass also increased in spring and summer throughout the water column, adding to organic enrichment. These increases were reflected by increased chl-a and phaeo-a in the lower column over the time series. The increases of large-celled phytoplankton and nanoflagellates were accompanied by elevated DON and bacteria, again in summer/autumn, indicating linkage between autotrophs and heterotrophic activity toward the end of the production season. The increasing biomasses during the time series of larger phytoplankton more prone to sedimentation suggested that deposition of biomass to sediments may have increased over the period, supported by the observations of increased sediment O 2 demand (Zeldis and Swaney 2018). Thus, enhanced organic matter deposition and enrichment of the sedimentary environment may have underlain decreased denitrification and increased availability of DIN to primary producers.
In the Firth, primary biomass and its mineralisation increased in the lower water column late in the production season, as shown by the elevated concentrations of chl-a, its degradation product (phaeo-a), and decreased O 2 in the lower water column between spring and summer (Fig. 9). These findings were consistent with the summer/autumn low O 2 regularly recorded in the lower half of the water column at the monitoring site over 1998-2013 (Zeldis and Swaney 2018). The results of Fig. 9 do not distinguish whether the deep biomass was the result of increased export from the upper water column or simply decreased production in surface waters and increased production at depth. However, while phytoplankton biomass has little seasonal variation in the Firth (Table 4), integrated net primary production (NPP) increases between spring and summer in the inner Hauraki Gulf and Firth (Bury et al. 2012;Gall and Zeldis 2011). This suggests that biomass turnover and export from the euphotic zone is likely to increase between spring and summer. This conforms with observations in Yangtze River estuary (China: Zhou et al. 2020), where low O 2 in bottom waters was driven by respiration of diatom-rich phytoplankton settling through the pycnocline in summer/autumn. In the Firth, the seasonal dynamics of biomass turnover is an important topic for future research, because of its association with expression of ecological stressors (acidification and hypoxia) as described by Zeldis et al. (2022). The increasing trends over 1998-2013 for larger phytoplankton taxa prone to sedimentation indicated potential for increased deposition of organic matter that could explain the biogeochemical changes (reduced denitrification and increased sediment O 2 demand) observed over the period (Zeldis and Swaney 2018).
Unlike the trends in DIN and microplankton, the 1998-2013 time series of water column O 2 at the monitoring site showed no important trend, with 'increasing trend about as likely as not' (not shown). The lack of consistent trend in water column O 2 could have arisen from the complex interacting factors that affect O 2 concentrations in the coastal sea, including stratification, air-sea exchange and degree of organic matter production (Hagy et al. 2004;Scully 2016), as was evident in the variable chl-a, stratification and O 2 we observed between autumn 2010 and autumn 2013 (Fig. 9). Irrespective, O 2 levels during 1998-2013 at the Firth site often reached 60% O 2 saturation and occasionally 40% (Zeldis and Swaney 2018). Thus, while the DIN and microplankton increases observed over 1998-2013 did not result in a worsening of hypoxia, they did indicate a furthering of the enrichment the Firth has received over multiple decades (Introduction) that drives frequent expression of low O 2 .

Summary and Conclusion
Our study has used a combination of time series and GAM analyses showing the key role of N enrichment in driving microplankton biomass increase over a decadal-scale time series, with weaker effects of stratification operating at shorter time scales. Important increases occurred for larger phytoplankton > 2 µm, including diatoms (dominated by increases of large centrics but also including increased rare but toxic Pseudo-nitzschia) and large increases of nanoflagellates. Dinoflagellates decreased, such that the community shifted from dinoflagellate to diatom/nanoflagellate-dominance. Among autotrophs, increased N:P may have favoured large diatoms and nanoflagellates with respect to dinoflagellates, while among heterotrophs, increases in bacteria and small ciliates were potentially favoured through N enrichment and food-web processes.
Because the enrichment was expressed not only as changed biomass but also changed community composition, there was potential for reduced nutritional quality of the phytoplankton for higher trophic levels, given the decrease of potentially more nutritious dinoflagellate biomass. The substantial increases of phytoplankton > 2 µm and large centric diatoms in particular indicated potential for enhanced sedimentation of organic matter, possibly explaining decreases in sediment denitrification recorded between early and late stages of the time series (Zeldis and Swaney 2018) and the source of the DIN increase over the time series.
Seasonality of the primary production cycle and nutrient limitation in the Firth underpin a transition from net-autotrophy in spring to early summer, toward net-heterotrophy in late summer and autumn (Zeldis et al. 2022) with late-season expression of depressed O 2 and pH linked to mineralisation of biomass. The observations made here of increased N and microplankton indicate a continuation of enrichment the Firth has received over decades and support previous conclusions that plankton at the base of the Firth food web are pivotal in the expression of ecosystem stressors (Zeldis and Swaney 2018). Our study thus shows the value of long-term records of microplankton population change, in identifying and explaining important drivers of coastal ecosystem health.
Acknowledgements We thank scientific staff and crews of NIWA R.V.s Kaharoa, Tangaroa and Western Work Boats Star Keys who enabled the ocean-going data, the NIWA Water Quality Laboratory for sample analysis, Jochen Bind for database administration and Scott Nodder for manuscript review. We thank journal reviewers and editors.
Funding Project support was from the New Zealand Ministry of Business, Innovation and Employment (MBIE)-funded 'Coasts and Oceans' Outcome-Based Investment contract C01X0501 and NIWA Strategic Science Investment Fund project 'Ocean Flows and Productivity' (COOF2202). Funds for open access fees were from 'Ocean Flows and Productivity' (COOF2202).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.