Phytoplankton as indicators of global warming?

Terrestrial plants are sensitive indicators of global warming because their annual cycles of growth and senescence are changing as warming proceeds. Single celled algae are distinct life forms capable of population bursts in any season, so there is uncertainty about phytoplankton phenology as a comparable indicator of global warming. We analyzed 4+ decades of monthly chlorophyll a measurements at two sites in San Francisco Bay and found abrupt shifts during summer months leading to a 48‐day advance in the annual pattern of chlorophyll‐a accumulation at one site and a 36‐day delay at the other. These large phenological changes were not associated with changing temperature, but they were associated with changes in top–down control by bivalve filter feeders as biological communities were restructured by (1) introduction of a non‐native clam, and (2) a shift in atmospheric forcing of the NE Pacific. This study illustrates that changes in phytoplankton phenology are not necessarily responses to or indicators of global warming. However, they can be indicators of human disturbances and natural climate oscillations having effects large enough to mask the effect of climate warming.

Terrestrial plants at mid and high latitudes have a canonical annual cycle of turning on (leafing out, flowering) in spring and turning off in autumn as responses to the annual cycles of temperature, precipitation, and solar radiation.These cycles are changing as the planet warms, and the most common change is a steady advancement of spring growth and extension of the growing season of terrestrial plants (Inouye 2022).These phenological changes provided early evidence that the Earth's biosphere is being modified by anthropogenic climate change (Parmesan and Yohe 2003).Changes in plant seasonal patterns have broad implications: they can reshape communities, alter the water cycle, increase annual primary production, facilitate establishment of non-native species, and have feedbacks to the climate system (Piao et al. 2019).Plants are sensitive indicators of changing climate because their life histories follow regular, repeatable annual cycles of growth, maturation, and senescence cued to the annual climate cycle.
Half of global primary production is phytoplankton photosynthesis (Field et al. 1998), and some observational records show phenological changes that parallel those of terrestrial plants-that is, monotonic trends of earlier seasonal blooms synchronous with local warming (e.g., Winder and Schindler 2004;Kahru and Elmgren 2014).However, the singlecelled algae are distinct life forms from plants because of their capacity to turn on in any season, fast growth rates, and tight biomass control by equally fast grazing losses (e.g., Hunter-Cevera et al. 2016).Unlike plants, there is no canonical annual cycle of phytoplankton biomass in lakes, estuaries, and oceans where blooms can occur any season and the annual cycles can change abruptly (Winder and Cloern 2010).Given this wide plasticity of annual cycles, there is uncertainty about changing phytoplankton phenology as an indicator of global warming.
Here, we analyzed monthly chlorophyll a (Chl a) time series from two regions of San Francisco Bay US to ask if phytoplankton at the land-sea interface follow the same phenology rules as terrestrial plants.In particular, we asked if 4+ decades of observations show: (1) advancement of the seasonal pattern, (2) synchrony of changing seasonal patterns across regions, and (3) changing seasonal patterns that can be attributed to warming.Our results suggest that, at least in some ecosystems, phytoplankton follow their own phenological rules.

Study site
San Francisco Bay is an estuarine ecosystem of the NE Pacific composed of two different estuary types (Fig. 1).The North Bay (NB) is a river-dominated estuary, and the South Bay (SB) is a marine-lagoon estuary influenced by its connection to the coastal NE Pacific.Early studies found distinct seasonal patterns of phytoplankton biomass, measured as chlorophyll-a concentration (hereafter, chl), in the two regions.The NB pattern was low chl in spring resulting from fast downstream transport during the high-inflow season, and peak chl in summer attributed to estuarine circulation that concentrates particles (including diatoms) in the landward NB during the low-inflow season (Cloern et al. 1983).The SB pattern was the opposite, with peak chl in spring and low chl in summer (Cloern 1982).The SB spring blooms developed North Bay during periods of salinity stratification, and the low summer chl was attributed to grazing control by bivalve filter feeders whose biomass and filtration rates increase during the warm season (Cloern 1982).From these studies, we learned that San Francisco Bay is highly enriched in nitrogen and phosphorus, nutrients are rarely limiting, and other factors (light limitation, grazing, and transport) drive phytoplankton variability (Cloern et al. 2020).

Data sources
San Francisco Bay has been a site of research and monitoring sustained over multiple decades by two programs: the Interagency Ecological Program (IEP) has conducted monthly sampling in NB since 1975 (Battey and Perry 2022), and the U.S. Geological Survey (USGS) began studies in 1969 (Conomos 1979) that include monthly or higher-frequency sampling in SB.From these programs, we compiled surface measurements of water quality variables made from shipbased sampling at three sites in NB by the IEP and at 12 sites in SB by USGS (Fig. 1).Measurements included chl as a proxy for phytoplankton biomass, specific conductance (NB) or salinity (SB), water temperature, light penetration as Secchi depth (NB) or light extinction coefficient (SB), dissolved inorganic nitrogen (DIN; nitrate + nitrite + ammonium), and phosphorus.The records include 1952 individual chl measurements in NB from 1975 to 2021 and 9681 chl measurements in SB from 1980 to 2021.We aggregated data to produce and analyze time series of monthly mean variables across the sampling sites in each region.
The data we compiled and supporting metadata are available online (Cloern 2023).

Analytical methods
North Bay Chl concentrations were measured in NB using USEPA Standard method 10,200 h.Temperature and specific conductance were measured in recent decades with a YSI EXO2 multiparameter sonde.Ammonium was measured in filtered samples with EPA method 350.1.Nitrate + nitrite was measured in recent decades with Standard Method 4500-NO3-F, and dissolved inorganic phosphorus (DIP) with EPA method 365.1.

South Bay
Chl measurements in SB included laboratory analyses of discrete water samples and in-vivo fluorometers calibrated with the discrete measurements.Temperature was measured in recent decades with Sea-Bird Electronics SBE-3 and SBE-3plus sensors.Salinity was measured with Sea-Bird Electronics SBE-4 or SBE-4C conductivity sensors.Light extinction coefficients were computed from vertical profiles of photosynthetically available radiation measured with a Li-Cor Biosciences LI-192 quantum sensor.Samples for dissolved inorganic nutrient analyses were filtered through polycarbonate 0.4-μm pore size membrane filters, and since 2014 filtrates were analyzed with a Thermo Scientific Aquakem 600 automated discrete analyzer using methods of Fishman and Friedman (1989) for nitrite and phosphate, the method of Patton and Kryskalla (2003) for nitrate, and the Solorzano (1969) method for ammonium.Schraga and Cloern (2017) provide detailed histories of methods changes over time and intercalibrations across methods.
More detailed sampling and analytical methods for both programs are provided in Supporting Information S1.

Biota
The IEP Environmental Monitoring Program has monitored macrobenthos communities in NB since 1975, identifying and counting invertebrates collected with Ponor dredges (Wells 2022).We used count data of the non-native clam Potamocorbula amurensis at site D7 (Fig. 1) as an indicator of benthic top-down control of phytoplankton biomass in this region of NB.The California Department of Fish and Wildlife Bay Study Program (Baxter et al. 1999) samples demersal fish, shrimp, and crab populations with otter trawls at 10 sites in SB that span the USGS sampling transect.We used this record as a proxy for bivalve-predator abundance in SB.

Data analyses
Data analyses were done in the R environment (R_Core_Team 2021), primarily with package wql (Jassby and Cloem 2015).We used functions eof to identify seasonal modes of chl variability in each region, seasonTrend to measure chl trends for individual months, the Pettitt test (function pett) to identify change points in annual series, and phenoPhase to measure changes in the phases of annual chl cycles.We used Spearman rank correlations to identify potential drivers of chl variability and generalized additive models (function gam, method REML) with package mgcv version 1.8.42 (Wood 2017) to compare effect sizes of those drivers.

Seasonal modes of chlorophyll variability
We used two methods to identify seasonal modes of chl variability from the 1975-2021 time series from NB and the 1980-2021 time series from SB. First, we used empirical orthogonal function (EOF) analysis applied to the year x month chl matrices.The first principal components explained 51% of the overall chl variability in NB and 36% of variability in the SB.The coefficients (eigenvalues) of the first EOF in both regions were small for the early months of the year and large for later months (Fig. S1).Second, we calculated chl trends (μg chl L À1 yr À1 ) by month for each region.This analysis showed small or nonsignificant trends for months January through April, but large and significant negative trends in NB and positive trends in SB for the summer months (right panels, Fig. S1).Based on these two analyses, we defined a spring mode (small EOF1 coefficients, weak or no trends) as chl averaged

Cloern et al. Phytoplankton as climate-change indicators
over the months January-April (NB) and February-April (SB), and a summer mode (large EOF1 coefficients, large chl trends) as mean chl from July-October (NB) and June-September (SB).The trends of summer chl did not follow monotonic patterns, but instead resulted from abrupt shifts (Fig. 2a).The first decade of observation in NB showed a pattern of high summer chl (median 11.6 μg L À1 ) and then a downshift after 1986 to a median of 1.9 μg L À1 .Conversely, the first two decades of observation in SB showed a consistent pattern of low summer chl (median 2.5 μg L À1 ) and then a discrete upshift after 2000 when the median increased to 5.5 μg L À1 .

Changing phytoplankton phenology
We used three phenological indicators to learn how these abrupt changes in summer chl altered the overall seasonal patterns of chl.First, we computed the ratio of summer to spring chl as an index of the relative importance of the two seasons when chl can be high.The Pettit test identified a NB change point in that ratio after 1986, when the median decreased from 4.36 to 0.97 (Fig. 2b).Although there is large annual variability of seasonal patterns (Fig. 2b), these differences between eras are statistically significant ( p = 8 Â 10 À3 ).Thus, the NB was transformed abruptly from a system where the seasonal pattern was dominated by summer blooms to one where summer blooms have disappeared and there is no difference between spring and summer chl.In SB, the summer: spring chl ratio changed (p = 7 Â 10 À5 ) after 2003, when the median increased from 0.32 to 0.95 (Fig. 2b).The SB was transformed from a system where chl was three times higher in spring than summer, to one where spring and summer chl are nearly equal.The net result of these shifts has been a homogenization of seasonal chl patterns in both regions.
Second, we used a common index of plankton phenology based on the monthly accumulation of chl each year (Greve et al. 2005).We used wql function phenoPhase to calculate the fulcrum (or "center of gravity") as the date each year (in months) when the cumulative chl reached half the total annual cumulative chl.Ecosystems dominated by spring (summer) blooms have small (large) fulcrum values.In NB, there was a significant (p = 1.6 Â 10 À3 ) decrease of the median fulcrum from month 7.06 to month 5.61 after 1986 (Fig. 2c).In SB, there was a significant (p = 5 Â 10 À4 ) increase of the median fulcrum from month 4.57 to month 5.74 after 2003 (Fig. 2c).
Third, we used cumulative sum functions (R function cumsum) of mean monthly chl time series to compare seasonal patterns of chl accumulation before and after the change points shown in Fig. 2b,c.These accumulation patterns changed in both regions, but only after spring (Fig. 2d), reflecting the absence of trends in the spring modes.The fulcrums of these cumsum functions indicate that the mean chl accumulation pattern in NB has advanced 48 d (from July 21 to June 3) as a result of the lost summer blooms.Conversely, the seasonal pattern in SB has receded 36 d (from April 25 to May 31) because of increased summer chl.These large abrupt changes in phytoplankton seasonal patterns are distinct from the more gradual changes in terrestrial-plant phenology at rates of about 2-9 d per decade (Vitasse et al. 2022).
The cumsum functions also show that the loss of summer blooms in NB reduced mean annual cumulative chl from 91 to 24 μg L À1 , and the summer chl upshift in SB increased the mean cumulative annual chl from 55 to 71 μg L À1 (Fig. 2d).Phytoplankton primary production scales linearly with chl in this estuary (Cole and Cloern 1987), so these changes in annual chl imply comparable changes in annual primary production.The rho is the Spearman rank correlation coefficient.Bold values are statistically significant using a p criterion of 0.05.n is the number of years where data were available for each potential driver of summer chl variability.

Processes of changing phytoplankton phenology
The phenological changes reported here result from changes in summer chl.We used correlations to identify potential drivers of summer chl variability and generalized additive models (GAMs) to measure the deviance of summer chl associated with individual drivers.We considered water temperature as a proxy for global warming, specific conductance (NB) or salinity (SB) as proxies for freshwater input and its effects on transport, stratification and mixing, Secchi depth (NB) or extinction coefficient (SB) as proxies for turbidity and light limitation of photosynthesis, DIN and DIP as pools of N and P that potentially limit phytoplankton growth rate.Our analyses showed no correlations between summer chl and summer water temperature, salinity, or light penetration in either region (Table 1).Therefore, any responses of summer chl to global warming or changes in hydrology or turbidity are masked by other processes.We found significant correlations between summer chl and both DIN and DIP, but those correlations were negative (Table 1) and reflect the drawdown of inorganic nutrients as phytoplankton biomass builds during blooms (Peterson et al. 1985).Over the course of our studies, the median summer DIN and DIP in NB (25.6 μM N, 2.7 μM P) and SB (28.6 μM N, 6.2 μM P) exceeded concentrations that limit phytoplankton growth.
We next considered the top-down process of phytoplankton consumption by bivalve filter feeders that can play strong roles in regulating phytoplankton biomass in shallow systems such as the littoral regions of lakes (Pothoven and Vanderploeg 2021), coastal bays (Beukema and Cadée 1996), and estuaries (Petersen et al. 2008).We focus on this grazer community because the loss of summer chl in NB (Fig. 2a) immediately followed the 1986 introduction and explosive Bivalve Predator Index Log Summer Chl chl in NB as a function of temperature and P. amurensis absence or presence; dots are jittered for clarity.The temperature effect was nonsignificant, the Potamocorbula effect was significant, and the model explained 82% of summer chl deviance.(c) An index of summer abundance of five species of demersal fish and crustaceans that prey on bivalve mollusks in SB.The red vertical line indicates a significant change point after 1999.(d) Partial residuals of a GAM of log-transformed summer chl in SB as a function of temperature and the bivalve-predator index above.The temperature effect was nonsignificant, the bivalve predator effect was significant, and the model explained 66% of chl deviance.
population growth (Fig. 3a) of the Asian clam P. amurensis (Alpine and Cloern 1992).Grazing by that introduced clam generally exceeds zooplankton grazing rates in NB (Kimmerer and Thompson 2014).Conversely, the chl upshift in SB (Fig. 2c) followed a decrease in bivalve abundance following population increases of demersal fish and crustaceans that prey on bivalves.That trophic cascade occurred after a 1999 shift in atmospheric forcing of the North Pacific, indexed as sign changes in both the Pacific Decadal Oscillation and North Pacific Gyre Oscillation (Cloern et al. 2010).
We found a significant negative correlation between summer chl and summer abundance of P. amurensis in NB (Table 1).We assessed the importance of P. amurensis and temperature as drivers of chl variability using a simple GAM: log(chl) $ s (temperature) + s (P.amurensis abundance).We treated P. amurensis abundance as a categorical variable (present or absent) because summer chl has been persistently low after its introduction, regardless of population size.This model showed no significant temperature effect (p = 0.19), a significant (p < 2 Â 10 À16 ) effect of P. amurensis, and accounted for 82% of summer chl deviance (Fig. 3b).
There is no comparable long-term record of bivalve abundance in SB, but there are abundance records of fish and crustacean species that prey on bivalves (Baxter et al. 1999).We found significant positive correlations between summer chl in SB and summer abundances of five common species of bivalve predators: English sole (Parophrys vetulus), plainfin midshipman (Porichthys notatus), blacktail bay shrimp (Crangon nigricauda), Stimpson coastal shrimp (Heptacarpus stimpsoni), and Pacific rock crab (Cancer antennarius).Abundances of these species vary by orders of magnitude, so we computed an index of bivalve-predator abundance as the mean of standardized series (mean zero, standard deviation 1) for these five species.This index had a significant shift from negative (below-mean predator abundance) to positive after 1999 (Fig. 3c), 1 yr before the upshift in summer chl (Fig. 2a).A GAM (log(chl) $ s (temperature) + s (predator index)) showed no association of summer chl with temperature (p = 0.82), a significant (p < 2 Â 10 À16 ) bivalvepredator effect, and explained 66% of summer chl deviance in SB (Fig. 3d).
These analyses show that the changes in summer chl and resulting seasonal patterns of phytoplankton biomass were not associated with temperature but were associated with changing abundance of bivalve filter feeders or their predators.Chl changes only in summer may result from seasonal variability of bivalve abundance and filtration rates that are both highest in summer (Cloern 1982).

Ecosystem implications of seasonal chlorophyll changes
An ecological function of phytoplankton is to convert dissolved forms of reactive elements into particulate forms that are consumed, metabolized, and regenerated to dissolved forms.The rates of these processes are ultimately determined by the rate of phytoplankton primary production, the largest source of autochthonous organic carbon to San Francisco Bay (Jassby et al. 1993).As chl increases (decreases), the rates of energy production and consumption, higher-level production, and element cycling accelerate (decelerate).We compiled examples from studies in San Francisco Bay to illustrate the systemic responses of an estuarine ecosystem to seasonal chl variability (Table 2).The rates of ecosystem functions-primary production, nutrient regeneration from sediments, pelagic and benthic respiration, and net ecosystem metabolism-all covary with chl.The abundance and growth rates of consumers-bacteria, protistans, polychaetes, crustacean zooplankton, and bivalves covary with chl.The concentrations of reactive elements (O 2 , C, N, P, Si, Cd, Ni, Zn, Pb, and Hg), C and N content of sediments, biochemical composition of consumer organisms, and assimilation of trace metals into food webs all covary with chl.Therefore, the changing phytoplankton phenology reported here (Fig. 2) implies comparable changes in the seasonal patterns of ecosystem metabolism, biological communities, and biogeochemical processes.

A global perspective
Decades of observation across the globe reveal three general features ("rules") of terrestrial-plant phenology: (1) climate change has altered the annual cycles of terrestrial plants, and the most common pattern is unidirectional advancement of spring growth (Piao et al. 2019); (2) phenological changes are synchronous across large (subcontinental) spatial scales (Shestakova et al. 2016); and (3) a primary control is temperature.Our study illustrates how phytoplankton do not necessarily follow any of these rules.We showed large changes in seasonal patterns within one ecosystem that were: (1) bidirectional, (2) asynchronous, (3) and not associated with warming.
To be clear, there are examples such as Lake Washington US (Winder and Schindler 2004), the Baltic Sea (Kahru and Elmgren 2014), and western Scheldt Estuary (Kromkamp and Engeland 2009) where trends of earlier phytoplankton blooms are synchronous with warming trends.And, Synechococcus blooms in the New England shelf develop earlier during warm springs (Hunter-Cevera et al. 2016).However, there are other places like San Francisco Bay where changes in phytoplankton biomass and phenology have been driven primarily by processes other than temperature-regulated growth rate, such as: nutrient enrichment of Lake Windermere UK (Thackeray et al. 2008); oligotrophication of Lake Mjøsa Norway (Moe et al. 2022); introduction of Dreissena mussels to Lake Michigan US (Pothoven and Vanderploeg 2021); invasion of the clam Mya arenaria to the Ringkøbing Fjord, Denmark (Petersen et al. 2008); and fishing harvest of top predators in the Baltic Sea (Casini et al. 2008).
In other places, phytoplankton seasonal patterns are driven by variability of large-scale atmospheric forcing.
Examples include changing bloom patterns in: Narraganset Bay US synchronous with fluctuations of the North Atlantic Oscillation (Borkman and Smayda 2009); the North Pacific Ocean tied to the Pacific Decadal Oscillation (Chiba et al. 2012); and the North Atlantic Ocean synchronous with the Atlantic Multidecadal Oscillation (Edwards et al. 2022).The North Atlantic case is another example of bidirectional patterns of change as diatom abundance has increased in the northern sectors and decreased in the southern sectors.
Our contemporary understanding of phytoplankton ecology is that seasonal patterns are highly dynamic, variability is driven by many processes that alter phytoplankton growth and grazing rates, responses to global warming are detectable in some ecosystems, but they can be masked in others by responses to other human disturbances and multidecadal oscillations of the climate system.Perhaps the singular distinguishing feature ("rule") of phytoplankton phenology is that its biomass can change rapidly and any time when the balance between production and consumption is altered.

Fig. 1 .
Fig. 1.(a) Map of San Francisco Bay showing locations of three IEP sampling sites in the landward region of NB and 12 USGS sampling sites in SB.(b) Time series of log-transformed surface Chl a concentrations (μg L À1 ) in NB from measurements made on 1390 dates from January 1975 to December 2021.(c) Time series of log-transformed surface Chl a concentrations in SB from measurements made on 903 dates from January 1980 to December 2021.

Fig. 2 .
Fig. 2. (a) Bar plots showing annual variability of summer chl (μg L À1 ) in North and South San Francisco Bay.Vertical red lines identify change points separating eras having statistically significant (NB p = 1 Â 10 À4 ; SB p = 4 Â 10 À7 ) differences in median summer chl.(b) Annual variability in the summer:spring chl ratio.That ratio decreased in NB from median 4.36 (standard deviation, SD = 2.69) to 0.97 (SD = 0.36) after 1986.The summer:spring ratio increased in SB from median 0.32 (SD = 0.15) to 0.95 (SD = 0.30) after 2003.(c) Annual variability in the seasonal pattern of chl variability measured as the fulcrumthe date in months when cumulative chl reached 50% of total annual cumulative chl.The median fulcrum decreased in NB from month 7.06 (SD = 0.97) to month 5.61 (SD = 0.67) after 1986.The median fulcrum increased in SB from month 4.57 (SD = 0.73) to month 5.74 (SD = 0.69) after 2003.(d) Mean patterns of chl accumulation (μg L À1 ) by month, before and after the shifts shown in (b,c).Filled circles show changes in the mean fulcrums after those shifts.

Fig. 3 .
Fig. 3. (a) Mean summer abundance (per m 2 ) of the clam P. amurensis introduced to NB in 1986.(b) Model predictions of log transformed summer

Table 1 .
Correlations between summer chl and potential drivers of its annual variability in North and South San Francisco Bay.

Table 2 .
Examples of ecosystem, organism, and biogeochemical processes or properties that covary (positive or negative) with Chl a concentrations in San Francisco Bay.